printf("failed to allocate memory for `StP' (%g MB) %d.\n", bytes / (1024.0 * 1024.0),sizeof(struct st_particle_data));
endrun(1);
}
bytes_tot += bytes;
if(ThisTask == 0)
printf("\nAllocated %g MByte for star properties storage. %d\n\n", bytes_tot / (1024.0 * 1024.0), sizeof(struct st_particle_data));
}
#endif
}
/*! This routine frees the memory for the particle storage. Note: We don't
* actually bother to call it in the code... When the program terminats,
* the memory will be automatically freed by the operating system.
*/
void free_memory(void)
{
#ifdef GAS_ACCRETION
if(N_gas_acc > 0)
free(SphAcc);
if(NumPart_acc > 0)
free(Acc);
#endif
#ifdef HOT_HALO
if(NumPart_HaloAcc > 0)
free(HaloAcc);
#endif
#ifdef STELLAR_PROP
if(All.MaxPartStars > 0)
free(StP);
#endif
if(All.MaxPartSph > 0)
free(SphP);
if(All.MaxPart > 0)
free(P);
#ifdef COOLING_FCT_FROM_HDF5
freeMemory();
#endif
}
diff --git a/src/allvars.c b/src/allvars.c
index 367ed3e..91e1611 100644
--- a/src/allvars.c
+++ b/src/allvars.c
@@ -1,578 +1,594 @@
/*! \file allvars.c
* \brief provides instances of all global variables.
*/
#include <stdio.h>
#include <gsl/gsl_rng.h>
#include "tags.h"
#include "allvars.h"
#ifdef CHIMIE
int FE=0;
int METALS=0;
#endif
int SetMinTimeStepForActives=0;
int ThisTask; /*!< the rank of the local processor */
int NTask; /*!< number of processors */
int PTask; /*!< smallest integer such that NTask <= 2^PTask */
int NumPart; /*!< number of particles on the LOCAL processor */
int N_gas; /*!< number of gas particles on the LOCAL processor */
#if defined(SFR) || defined(STELLAR_PROP)
int N_stars; /*!< number of stars particle on the LOCAL processor */
#endif
#ifdef MULTIPHASE
int N_sph;
int N_sticky;
int N_stickyflaged;
int N_dark;
int NumColPotLocal; /*!< local number of potentially collisional particles */
int NumColPot; /*!< total number of potentially collisional particles */
int NumColLocal; /*!< local number of collisions */
int NumCol; /*!< total number of collisions */
int NumNoColLocal;
int NumNoCol;
#endif
#ifdef GAS_ACCRETION
int NumPart_acc;
int N_gas_acc;
#ifdef STELLAR_PROP
int N_stars_acc;
#endif
#endif
#ifdef HOT_HALO
int NumPart_HaloAcc;
#endif
long long Ntype[6]; /*!< total number of particles of each type */
int NtypeLocal[6]; /*!< local number of particles of each type */
int NumForceUpdate; /*!< number of active particles on local processor in current timestep */
int NumSphUpdate; /*!< number of active SPH particles on local processor in current timestep */
#ifdef CHIMIE
int NumStUpdate;
#endif
#ifdef TESSEL
int NumPTUpdate;
#endif
double CPUThisRun; /*!< Sums the CPU time for the process (current submission only) */
#ifdef SPLIT_DOMAIN_USING_TIME
double CPU_Gravity;
#endif
int RestartFlag; /*!< taken from command line used to start code. 0 is normal start-up from
initial conditions, 1 is resuming a run from a set of restart files, while 2
marks a restart from a snapshot file. */
char *Exportflag; /*!< Buffer used for flagging whether a particle needs to be exported to another process */
int *Ngblist; /*!< Buffer to hold indices of neighbours retrieved by the neighbour search routines */
int TreeReconstructFlag; /*!< Signals that a new tree needs to be constructed */
#ifdef SFR
int RearrangeParticlesFlag;/*!< Signals that particles must be rearanged */
#endif
int Flag_FullStep; /*!< This flag signals that the current step involves all particles */
gsl_rng *random_generator; /*!< the employed random number generator of the GSL library */
+#ifdef SIDM_ON
+gsl_rng *random_generator_SIDM; /*!< the employed random number generator of the GSL library here for SIDM calculations*/
+#endif
+
+
double RndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#ifdef SFR
double StarFormationRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef FEEDBACK_WIND
double FeedbackWindRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE
double ChimieRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
double ChimieKineticFeedbackRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef GAS_ACCRETION
double gasAccretionRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef HOT_HALO
double HaloRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef AB_TURB
//Ornstein-Uhlenbeck variables
double StOUVar;
double* StOUPhases;
gsl_rng* StRng;
//forcing field in fourie space
double* StAmpl;
double* StAka; //phases (real part)
double* StAkb; //phases (imag part)
double* StMode;
int StNModes;
//integertime StTPrev; (yr : ask ?)
int StTPrev;
double StSolWeightNorm;
#endif
#ifdef PY_INTERFACE
int NumPartQ;
int N_gasQ;
long long NtypeQ[6]; /*!< total number of particles of each type */
int NtypeLocalQ[6]; /*!< local number of particles of each type */
double DomainCornerQ[3]; /*!< gives the lower left corner of simulation volume */
double DomainCenterQ[3]; /*!< gives the center of simulation volume */
double DomainLenQ; /*!< gives the (maximum) side-length of simulation volume */
double DomainFacQ; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
int DomainMyStartQ; /*!< first domain mesh cell that resides on the local processor */
int DomainMyLastQ; /*!< last domain mesh cell that resides on the local processor */
int *DomainStartListQ; /*!< a table that lists the first domain mesh cell for all processors */
int *DomainEndListQ; /*!< a table that lists the last domain mesh cell for all processors */
double *DomainWorkQ; /*!< a table that gives the total "work" due to the particles stored by each processor */
int *DomainCountQ; /*!< a table that gives the total number of particles held by each processor */
int *DomainCountSphQ; /*!< a table that gives the total number of SPH particles held by each processor */
int *DomainTaskQ; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
peanokey *DomainKeyBufQ; /*!< this points to a buffer used during the exchange of particle data */
int NTopnodesQ; /*!< total number of nodes in top-level tree */
int NTopleavesQ; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
void *CommBufferQ; /*!< points to communication buffer, which is used in the domain decomposition, the
parallel tree-force computation, the SPH routines, etc. */
int NTopnodesQ; /*!< total number of nodes in top-level tree */
int NTopleavesQ; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
#endif
double DomainCorner[3]; /*!< gives the lower left corner of simulation volume */
double DomainCenter[3]; /*!< gives the center of simulation volume */
double DomainLen; /*!< gives the (maximum) side-length of simulation volume */
double DomainFac; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
int DomainMyStart; /*!< first domain mesh cell that resides on the local processor */
int DomainMyLast; /*!< last domain mesh cell that resides on the local processor */
int *DomainStartList; /*!< a table that lists the first domain mesh cell for all processors */
int *DomainEndList; /*!< a table that lists the last domain mesh cell for all processors */
double *DomainWork; /*!< a table that gives the total "work" due to the particles stored by each processor */
int *DomainCount; /*!< a table that gives the total number of particles held by each processor */
int *DomainCountSph; /*!< a table that gives the total number of SPH particles held by each processor */
int *DomainTask; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
int *DomainNodeIndex; /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */
FLOAT *DomainTreeNodeLen; /*!< this table gives for each leaf of the top-level tree the side-length of the corresponding node of the gravitational tree */
FLOAT *DomainHmax; /*!< this table gives for each leaf of the top-level tree the maximum SPH smoothing length among the particles of the corresponding node of the gravitational tree */
struct DomainNODE
*DomainMoment; /*!< this table stores for each node of the top-level tree corresponding node data from the gravitational tree */
peanokey *DomainKeyBuf; /*!< this points to a buffer used during the exchange of particle data */
peanokey *Key; /*!< a table used for storing Peano-Hilbert keys for particles */
peanokey *KeySorted; /*!< holds a sorted table of Peano-Hilbert keys for all particles, used to construct top-level tree */
int NTopnodes; /*!< total number of nodes in top-level tree */
int NTopleaves; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
struct topnode_data
*TopNodes; /*!< points to the root node of the top-level tree */
#ifdef PY_INTERFACE
struct topnode_data
*TopNodesQ; /*!< points to the root node of the top-level tree */
#endif
double TimeOfLastTreeConstruction; /*!< holds what it says, only used in connection with FORCETEST */
/* variables for input/output, usually only used on process 0 */
char ParameterFile[MAXLEN_FILENAME]; /*!< file name of parameterfile used for starting the simulation */
FILE *FdInfo; /*!< file handle for info.txt log-file. */
FILE *FdLog ; /*!< file handle for log.txt log-file. */
FILE *FdEnergy; /*!< file handle for energy.txt log-file. */
#ifdef SYSTEMSTATISTICS
FILE *FdSystem;
#endif
FILE *FdTimings; /*!< file handle for timings.txt log-file. */
FILE *FdCPU; /*!< file handle for cpu.txt log-file. */
+#ifdef SIDM_ON
+FILE *FdSIDMScatter; /*!< file handle for scatter.txt log-files (Different one for each core). */
+#endif
#ifdef FORCETEST
FILE *FdForceTest; /*!< file handle for forcetest.txt log-file. */
#endif
#ifdef SFR
FILE *FdSfr; /*!< file handle for sfr.txt log-file. */
#endif
#ifdef CHIMIE
FILE *FdChimie; /*!< file handle for chimie log-file. */
#ifdef CHIMIE_STATS
FILE *FdChimieStatsSNs; /*!< file handle for chimie stats-file. */
FILE *FdChimieStatsGas; /*!< file handle for chimie stats-file. */
#endif
#ifdef FOF
FILE *FdFOF_Chimie; /*!< file handle for fof chimie log-file. */
#endif
#endif
#ifdef MULTIPHASE
FILE *FdPhase; /*!< file handle for phase.txt log-file. */
FILE *FdSticky; /*!< file handle for sticky.txt log-file. */
#endif
#ifdef AGN_ACCRETION
FILE *FdAccretion; /*!< file handle for accretion.txt log-file. */
#endif
#ifdef BONDI_ACCRETION
FILE *FdBondi; /*!< file handle for bondi.txt log-file. */
#endif
#ifdef BUBBLES
FILE *FdBubble; /*!< file handle for bubble.txt log-file. */
#endif
#ifdef GAS_ACCRETION
FILE *FdGasAccretion; /*!< file handle for gas_accretion.txt log-file. */
#endif
#ifdef PERIODICOUTER
FILE *FdTrace; /*!< file handle for trace.txt log-file. */
#endif
double DriftTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological drift factors */
double GravKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for gravitational forces */
double HydroKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for hydrodynmical forces */
#ifdef COSMICTIME
double CosmicTimeTable[COSMICTIME_TABLE_LENGTH]; /*!< table for the computation of cosmic time */
double FullCosmicTimeTable[COSMICTIME_TABLE_LENGTH]; /*!< table for the computation of cosmic time */
double FullCosmicTimeTableInv[COSMICTIME_TABLE_LENGTH]; /*!< table for the computation of cosmic time */
#endif
void *CommBuffer; /*!< points to communication buffer, which is used in the domain decomposition, the
parallel tree-force computation, the SPH routines, etc. */
/*! This structure contains data which is the SAME for all tasks (mostly code parameters read from the
* parameter file). Holding this data in a structure is convenient for writing/reading the restart file, and
* it allows the introduction of new global variables in a simple way. The only thing to do is to introduce
* them into this structure.
*/
struct global_data_all_processes
All;
/*! This structure holds all the information that is
* stored for each particle of the simulation.
*/
struct particle_data
*P, /*!< holds particle data on local processor */
*DomainPartBuf; /*!< buffer for particle data used in domain decomposition */
#ifdef PY_INTERFACE
struct particle_data
*Q, /*!< holds particle data on local processor */
*DomainPartBufQ; /*!< buffer for particle data used in domain decomposition */
#endif
/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
* variables.
*/
struct sph_particle_data
*SphP, /*!< holds SPH particle data on local processor */
*DomainSphBuf; /*!< buffer for SPH particle data in domain decomposition */
#ifdef FOF
struct fof_particles_in
*FoF_P,*FoF_P_local;
struct fof_groups_data
*Groups,*GroupsGet,*GroupsRecv;
int *Tot_HeadID_list;
int Tot_Ngroups;
int Tot_Ntgroups;
int Tot_Npgroups;
int Tot_Ntsfrgroups;
int Ngroups; /* true and pseudo groups */
int Ntgroups; /* true groups */
int Npgroups; /* pseudo groups */
int Ntsfrgroups; /* true groups that forms stars */
int Nsfrgroups; /* local number of groups flagged for sfr */
int Tot_Nsfrgroups; /* total number of groups flagged for sfr */
#endif
#ifdef GAS_ACCRETION
struct acc_particle_data
*Acc;
struct gas_acc_particle_data
*SphAcc;
#endif
#ifdef HOT_HALO
struct haloacc_pos_data
*HaloAcc;
#endif
#ifdef TRACE_ACC
struct pot_index
*PotInd;
#endif
#ifdef TRACE_ACC
struct pot_index
*PotInd;
#endif
#ifdef PY_INTERFACE
struct sph_particle_data
*SphQ, /*!< holds SPH particle data on local processor */
*DomainSphBufQ; /*!< buffer for SPH particle data in domain decomposition */
#endif
#ifdef STELLAR_PROP
/* the following struture holds data that is stored for each SPH particle in addition to the collisionless
* variables.
*/
struct st_particle_data
*StP, /*!< holds ST particle data on local processor */
*DomainStBuf; /*!< buffer for ST particle data in domain decomposition */
#endif
/* Variables for Tree
*/
int MaxNodes; /*!< maximum allowed number of internal nodes */
int Numnodestree; /*!< number of (internal) nodes in each tree */
struct NODE
*Nodes_base, /*!< points to the actual memory allocted for the nodes */
*Nodes; /*!< this is a pointer used to access the nodes which is shifted such that Nodes[All.MaxPart]
gives the first allocated node */
int *Nextnode; /*!< gives next node in tree walk */
int *Father; /*!< gives parent node in tree */
struct extNODE /*!< this structure holds additional tree-node information which is not needed in the actual gravity computation */
*Extnodes_base, /*!< points to the actual memory allocted for the extended node information */
*Extnodes; /*!< provides shifted access to extended node information, parallel to Nodes/Nodes_base */
/*! Header for the standard file format.
*/
struct io_header
header; /*!< holds header for snapshot files */
#ifdef CHIMIE_EXTRAHEADER
/*! Header for the chimie part.
*/
struct io_chimie_extraheader
chimie_extraheader;
#endif
char Tab_IO_Labels[IO_NBLOCKS][4]; /*<! This table holds four-byte character tags used for fileformat 2 */
/* global state of system, used for global statistics
*/
struct state_of_system
SysState; /*<! Structure for storing some global statistics about the simulation. */
/*! This structure contains data related to the energy budget.
These values are different for each task. It need to be stored
in the restart flag.
*/
struct local_state_of_system
LocalSysState; /*<! Structure for storing some local statistics about the simulation. */
/* Various structures for communication
*/
struct gravdata_in
*GravDataIn, /*!< holds particle data to be exported to other processors */
*GravDataGet, /*!< holds particle data imported from other processors */
*GravDataResult, /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
*GravDataOut; /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */
struct gravdata_index
*GravDataIndexTable; /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */
struct densdata_in
*DensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensDataGet; /*!< holds imported particle data for SPH density computation */
struct densdata_out
*DensDataResult, /*!< stores the locally computed SPH density results for imported particles */
*DensDataPartialResult; /*!< imported partial SPH density results from other processors */
struct hydrodata_in
*HydroDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*HydroDataGet; /*!< holds imported particle data for SPH hydro-force computation */
struct hydrodata_out
*HydroDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*HydroDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
#ifdef MULTIPHASE
struct stickydata_in
*StickyDataIn,
*StickyDataGet;
struct stickydata_out
*StickyDataResult,
*StickyDataPartialResult;
struct Sticky_index *StickyIndex;
#endif
#ifdef CHIMIE
struct chimiedata_in
*ChimieDataIn, /*!< holds particle data for Chimie computation to be exported to other processors */
*ChimieDataGet; /*!< holds imported particle data for Chimie computation */
struct chimiedata_out
*ChimieDataResult, /*!< stores the locally computed Chimie results for imported particles */
*ChimieDataPartialResult; /*!< imported partial Chimie results from other processors */
struct starsdensdata_in
*StarsDensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*StarsDensDataGet; /*!< holds imported particle data for SPH density computation */
struct starsdensdata_out
*StarsDensDataResult, /*!< stores the locally computed SPH density results for imported particles */
*StarsDensDataPartialResult; /*!< imported partial SPH density results from other processors */
#endif
#ifdef DISSIPATION_FORCES
struct dissipationforcesdata_in
*DissipationForcesDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*DissipationForcesDataGet; /*!< holds imported particle data for SPH hydro-force computation */
struct dissipationforcesdata_out
*DissipationForcesDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*DissipationForcesDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
#endif /* DISSIPATION_FORCES */
#ifdef FOF
struct FOFdata_in
*FOFDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*FOFDataGet; /*!< holds imported particle data for SPH hydro-force computation */
struct FOFdata_out
*FOFDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*FOFDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
struct FOFgravdata_in
*FOFGravDataIn, /*!< holds particle data to be exported to other processors */
*FOFGravDataGet, /*!< holds particle data imported from other processors */
*FOFGravDataResult, /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
*FOFGravDataOut; /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */
struct FOFgravdata_index
*FOFGravDataIndexTable;
#endif /* FOF */
#ifdef TESSEL
struct ghostdata_in
*GhostDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*GhostDataGet; /*!< holds imported particle data for SPH density computation */
struct ghostdata_out
*GhostDataResult, /*!< stores the locally computed SPH density results for imported particles */
*GhostDataPartialResult; /*!< imported partial SPH density results from other processors */
//struct ghost_particle_data
// *gP;
int NumgPart;
#endif /* TESSEL */
#ifdef SYNCHRONIZE_NGB_TIMESTEP
struct SynchroinzeNgbTimestepdata_in
*SynchroinzeNgbTimestepDataIn,
*SynchroinzeNgbTimestepDataGet;
#endif
#ifdef PY_INTERFACE
struct sphdata_in
*SphDataIn,
*SphDataGet;
struct sphdata_out
*SphDataResult,
*SphDataPartialResult;
struct denssphdata_in
*DensSphDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensSphDataGet; /*!< holds imported particle data for SPH density computation */
struct denssphdata_out
*DensSphDataResult, /*!< stores the locally computed SPH density results for imported particles */
*DensSphDataPartialResult; /*!< imported partial SPH density results from other processors */
#endif
+#ifdef SIDM_ON
+struct SIDMdata_in
+*SIDMDataIn, /*!< holds particle data for SIDM scattering to be exported to other processors */
+*SIDMDataGet; /*!< holds imported particle data for SIDM scattering */
+struct SIDMdata_out
+*SIDMDataResult, /*!< stores the locally computed SPH density results for imported particles */
+*SIDMDataPartialResult; /*!< imported partial SPH density results from other processors */
+#endif
diff --git a/src/allvars.h b/src/allvars.h
index f450339..0991118 100644
--- a/src/allvars.h
+++ b/src/allvars.h
@@ -1,2668 +1,2717 @@
/*! \file allvars.h
* \brief declares global variables.
*
* This file declares all global variables. Further variables should be added here, and declared as
* 'extern'. The actual existence of these variables is provided by the file 'allvars.c'. To produce
* 'allvars.c' from 'allvars.h', do the following:
*
* - Erase all #define's, typedef's, and enum's
* - add #include "allvars.h", delete the #ifndef ALLVARS_H conditional
* - delete all keywords 'extern'
* - delete all struct definitions enclosed in {...}, e.g.
#define GADGETVERSION "2.0" /*!< code version string */
#define TIMEBASE (1<<28) /*!< The simulated timespan is mapped onto the integer interval [0,TIMESPAN],
* where TIMESPAN needs to be a power of 2. Note that (1<<28) corresponds to 2^29
*/
#define MAXTOPNODES 200000 /*!< Maximum number of nodes in the top-level tree used for domain decomposition */
typedef long long peanokey; /*!< defines the variable type used for Peano-Hilbert keys */
#define BITS_PER_DIMENSION 18 /*!< Bits per dimension available for Peano-Hilbert order.
Note: If peanokey is defined as type int, the allowed maximum is 10.
If 64-bit integers are used, the maximum is 21 */
#define PEANOCELLS (((peanokey)1)<<(3*BITS_PER_DIMENSION)) /*!< The number of different Peano-Hilbert cells */
#define RNDTABLE 3000 /*!< gives the length of a table with random numbers, refreshed at every timestep.
This is used to allow application of random numbers to a specific particle
in a way that is independent of the number of processors used. */
#define MAX_REAL_NUMBER 1e37
#define MIN_REAL_NUMBER 1e-37
#define MAXLEN_FILENAME 100 /*!< Maximum number of characters for filenames (including the full path) */
#ifdef ISOTHERM_EQS
#define GAMMA (1.0) /*!< index for isothermal gas */
#else
#define GAMMA (5.0/3) /*!< adiabatic index of simulated gas */
#endif
#define GAMMA_MINUS1 (GAMMA-1)
#define HYDROGEN_MASSFRAC 0.76 /*!< mass fraction of hydrogen, relevant only for radiative cooling */
/* Some physical constants in cgs units */
#define GRAVITY 6.672e-8 /*!< Gravitational constant (in cgs units) */
#define SOLAR_MASS 1.989e33
#define SOLAR_LUM 3.826e33
#define RAD_CONST 7.565e-15
#define AVOGADRO 6.0222e23
#define BOLTZMANN 1.3806e-16
#define GAS_CONST 8.31425e7
#define C 2.9979e10
#define PLANCK 6.6262e-27
#define CM_PER_MPC 3.085678e24
#define PROTONMASS 1.6726e-24
#define ELECTRONMASS 9.10953e-28
#define THOMPSON 6.65245e-25
#define ELECTRONCHARGE 4.8032e-10
#define HUBBLE 3.2407789e-18 /* in h/sec */
#define YEAR_IN_SECOND 31536000.0 /* year in sec */
#define FEH_SOLAR 0.00181 /* used only if cooling with metal is on and chimie is off */
#define KPC_IN_CM 3.085678e+21
#define PI 3.1415926535897931
#define TWOPI 6.2831853071795862
/* Some conversion factors */
#define SEC_PER_MEGAYEAR 3.155e13
#define SEC_PER_YEAR 3.155e7
#ifndef ASMTH
#define ASMTH 1.25 /*!< ASMTH gives the scale of the short-range/long-range force split in units of FFT-mesh cells */
#endif
#ifndef RCUT
#define RCUT 4.5 /*!< RCUT gives the maximum distance (in units of the scale used for the force split) out to
which short-range forces are evaluated in the short-range tree walk. */
#endif
#define MAX_NGB 20000 /*!< defines maximum length of neighbour list */
#define MAXLEN_OUTPUTLIST 10000 /*!< maxmimum number of entries in list of snapshot output times */
#define DRIFT_TABLE_LENGTH 1000 /*!< length of the lookup table used to hold the drift and kick factors */
#ifdef COSMICTIME
#define COSMICTIME_TABLE_LENGTH 1000 /*!< length of the lookup table used for the cosmic time computation */
#endif
#define MAXITER 1000 /*!< maxmimum number of steps for SPH neighbour iteration */
#ifdef DOUBLEPRECISION /*!< If defined, the variable type FLOAT is set to "double", otherwise to FLOAT */
#define FLOAT double
#else
#define FLOAT float
#endif
#ifndef TWODIMS
#define NUMDIMS 3 /*!< For 3D-normalized kernel */
#define KERNEL_COEFF_1 2.546479089470 /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 15.278874536822
#define KERNEL_COEFF_3 45.836623610466
#define KERNEL_COEFF_4 30.557749073644
#define KERNEL_COEFF_5 5.092958178941
#define KERNEL_COEFF_6 (-15.278874536822)
#define NORM_COEFF 4.188790204786 /*!< Coefficient for kernel normalization. Note: 4.0/3 * PI = 4.188790204786 */
#else
#define NUMDIMS 2 /*!< For 2D-normalized kernel */
#define KERNEL_COEFF_1 (5.0/7*2.546479089470) /*!< Coefficients for SPH spline kernel and its derivative */
#define KERNEL_COEFF_2 (5.0/7*15.278874536822)
#define KERNEL_COEFF_3 (5.0/7*45.836623610466)
#define KERNEL_COEFF_4 (5.0/7*30.557749073644)
#define KERNEL_COEFF_5 (5.0/7*5.092958178941)
#define KERNEL_COEFF_6 (5.0/7*(-15.278874536822))
#define NORM_COEFF M_PI /*!< Coefficient for kernel normalization. */
#endif
#ifdef MULTIPHASE
#define GAS_SPH 0
#define GAS_STICKY 1
#define GAS_DARK 2
#endif
#if defined(SFR) || defined(STELLAR_PROP)
#define ST 1
#endif
#ifdef CHIMIE
#define NELEMENTS 10
#define MAXNELEMENTS 64
#endif
#ifdef COOLING
#define COOLING_NMETALICITIES 9
#define COOLING_NTEMPERATURES 171
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
#define VELOCITY_DISPERSION_SIZE 3
#endif
#ifdef TRACE_ACC
#define NUMTRACEAVG 64
#endif
#ifdef CHIMIE
extern int FE;
extern int METALS;
#endif
extern int SetMinTimeStepForActives;
extern int ThisTask; /*!< the rank of the local processor */
extern int NTask; /*!< number of processors */
extern int PTask; /*!< smallest integer such that NTask <= 2^PTask */
extern int NumPart; /*!< number of particles on the LOCAL processor */
extern int N_gas; /*!< number of gas particles on the LOCAL processor */
#if defined(SFR) || defined(STELLAR_PROP)
extern int N_stars; /*!< number of stars particle on the LOCAL processor */
#endif
#ifdef MULTIPHASE
extern int N_sph;
extern int N_sticky;
extern int N_stickyflaged;
extern int N_dark;
extern int NumColPotLocal; /*!< local number of potentially collisional particles */
extern int NumColPot; /*!< total number of potentially collisional particles */
extern int NumColLocal; /*!< local number of collisions */
extern int NumCol; /*!< total number of collisions */
extern int NumNoColLocal;
extern int NumNoCol;
#endif
#ifdef GAS_ACCRETION
extern int NumPart_acc;
extern int N_gas_acc;
#ifdef STELLAR_PROP
extern int N_stars_acc;
#endif
#endif
#ifdef HOT_HALO
extern int NumPart_HaloAcc;
#endif
extern long long Ntype[6]; /*!< total number of particles of each type */
extern int NtypeLocal[6]; /*!< local number of particles of each type */
extern int NumForceUpdate; /*!< number of active particles on local processor in current timestep */
extern int NumSphUpdate; /*!< number of active SPH particles on local processor in current timestep */
#ifdef CHIMIE
extern int NumStUpdate;
#endif
#ifdef TESSEL
extern int NumPTUpdate;
#endif
extern double CPUThisRun; /*!< Sums the CPU time for the process (current submission only) */
#ifdef SPLIT_DOMAIN_USING_TIME
extern double CPU_Gravity;
#endif
extern int RestartFlag; /*!< taken from command line used to start code. 0 is normal start-up from
initial conditions, 1 is resuming a run from a set of restart files, while 2
marks a restart from a snapshot file. */
extern char *Exportflag; /*!< Buffer used for flagging whether a particle needs to be exported to another process */
extern int *Ngblist; /*!< Buffer to hold indices of neighbours retrieved by the neighbour search routines */
extern int TreeReconstructFlag; /*!< Signals that a new tree needs to be constructed */
#ifdef SFR
extern int RearrangeParticlesFlag;/*!< Signals that particles must be rearanged */
#endif
extern int Flag_FullStep; /*!< This flag signals that the current step involves all particles */
extern gsl_rng *random_generator; /*!< the employed random number generator of the GSL library */
+#ifdef SIDM_ON
+extern gsl_rng *random_generator_SIDM; /*!< the employed random number generator of the GSL library here for SIDM calculations*/
+#endif
+
extern double RndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#ifdef SFR
extern double StarFormationRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef FEEDBACK_WIND
extern double FeedbackWindRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE
extern double ChimieRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
extern double ChimieKineticFeedbackRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef GAS_ACCRETION
extern double gasAccretionRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef HOT_HALO
extern double HaloRndTable[RNDTABLE]; /*!< Hold a table with random numbers, refreshed every timestep */
#endif
#ifdef AB_TURB
//Ornstein-Uhlenbeck variables
extern double StOUVar;
extern double* StOUPhases;
extern gsl_rng* StRng;
//forcing field in fourie space
extern double* StAmpl;
extern double* StAka; //phases (real part)
extern double* StAkb; //phases (imag part)
extern double* StMode;
extern int StNModes;
//integertime StTPrev; (yr : ask ?)
extern int StTPrev;
extern double StSolWeightNorm;
#endif
#ifdef PY_INTERFACE
extern int NumPartQ;
extern int N_gasQ;
extern long long NtypeQ[6]; /*!< total number of particles of each type */
extern int NtypeLocalQ[6]; /*!< local number of particles of each type */
extern double DomainCornerQ[3]; /*!< gives the lower left corner of simulation volume */
extern double DomainCenterQ[3]; /*!< gives the center of simulation volume */
extern double DomainLenQ; /*!< gives the (maximum) side-length of simulation volume */
extern double DomainFacQ; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
extern int DomainMyStartQ; /*!< first domain mesh cell that resides on the local processor */
extern int DomainMyLastQ; /*!< last domain mesh cell that resides on the local processor */
extern int *DomainStartListQ; /*!< a table that lists the first domain mesh cell for all processors */
extern int *DomainEndListQ; /*!< a table that lists the last domain mesh cell for all processors */
extern double *DomainWorkQ; /*!< a table that gives the total "work" due to the particles stored by each processor */
extern int *DomainCountQ; /*!< a table that gives the total number of particles held by each processor */
extern int *DomainCountSphQ; /*!< a table that gives the total number of SPH particles held by each processor */
extern int *DomainTaskQ; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
extern peanokey *DomainKeyBufQ; /*!< this points to a buffer used during the exchange of particle data */
extern int NTopnodesQ; /*!< total number of nodes in top-level tree */
extern int NTopleavesQ; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
extern void *CommBufferQ; /*!< points to communication buffer, which is used in the domain decomposition, the
parallel tree-force computation, the SPH routines, etc. */
#endif
extern double DomainCorner[3]; /*!< gives the lower left corner of simulation volume */
extern double DomainCenter[3]; /*!< gives the center of simulation volume */
extern double DomainLen; /*!< gives the (maximum) side-length of simulation volume */
extern double DomainFac; /*!< factor used for converting particle coordinates to a Peano-Hilbert mesh covering the simulation volume */
extern int DomainMyStart; /*!< first domain mesh cell that resides on the local processor */
extern int DomainMyLast; /*!< last domain mesh cell that resides on the local processor */
extern int *DomainStartList; /*!< a table that lists the first domain mesh cell for all processors */
extern int *DomainEndList; /*!< a table that lists the last domain mesh cell for all processors */
extern double *DomainWork; /*!< a table that gives the total "work" due to the particles stored by each processor */
extern int *DomainCount; /*!< a table that gives the total number of particles held by each processor */
extern int *DomainCountSph; /*!< a table that gives the total number of SPH particles held by each processor */
extern int *DomainTask; /*!< this table gives for each leaf of the top-level tree the processor it was assigned to */
extern int *DomainNodeIndex; /*!< this table gives for each leaf of the top-level tree the corresponding node of the gravitational tree */
extern FLOAT *DomainTreeNodeLen; /*!< this table gives for each leaf of the top-level tree the side-length of the corresponding node of the gravitational tree */
extern FLOAT *DomainHmax; /*!< this table gives for each leaf of the top-level tree the maximum SPH smoothing length among the particles of the corresponding node of the gravitational tree */
extern struct DomainNODE
{
FLOAT s[3]; /*!< center-of-mass coordinates */
FLOAT vs[3]; /*!< center-of-mass velocities */
FLOAT mass; /*!< mass of node */
#ifdef STELLAR_FLUX
FLOAT starlum; /*!< star luminosity of node */
#endif
#ifdef UNEQUALSOFTENINGS
#ifndef ADAPTIVE_GRAVSOFT_FORGAS
int bitflags; /*!< this bit-field encodes the particle type with the largest softening among the particles of the nodes, and whether there are particles with different softening in the node */
#else
FLOAT maxsoft; /*!< hold the maximum gravitational softening of particles in the
node if the ADAPTIVE_GRAVSOFT_FORGAS option is selected */
#endif
#endif
}
*DomainMoment; /*!< this table stores for each node of the top-level tree corresponding node data from the gravitational tree */
extern peanokey *DomainKeyBuf; /*!< this points to a buffer used during the exchange of particle data */
extern peanokey *Key; /*!< a table used for storing Peano-Hilbert keys for particles */
extern peanokey *KeySorted; /*!< holds a sorted table of Peano-Hilbert keys for all particles, used to construct top-level tree */
extern int NTopnodes; /*!< total number of nodes in top-level tree */
extern int NTopleaves; /*!< number of leaves in top-level tree. Each leaf can be assigned to a different processor */
extern struct topnode_data
{
int Daughter; /*!< index of first daughter cell (out of 8) of top-level node */
int Pstart; /*!< for the present top-level node, this gives the index of the first node in the concatenated list of topnodes collected from all processors */
int Blocks; /*!< for the present top-level node, this gives the number of corresponding nodes in the concatenated list of topnodes collected from all processors */
int Leaf; /*!< if the node is a leaf, this gives its number when all leaves are traversed in Peano-Hilbert order */
peanokey Size; /*!< number of Peano-Hilbert mesh-cells represented by top-level node */
peanokey StartKey; /*!< first Peano-Hilbert key in top-level node */
long long Count; /*!< counts the number of particles in this top-level node */
}
#ifdef PY_INTERFACE
*TopNodesQ,
#endif
*TopNodes; /*!< points to the root node of the top-level tree */
extern double TimeOfLastTreeConstruction; /*!< holds what it says, only used in connection with FORCETEST */
/* variables for input/output, usually only used on process 0 */
extern char ParameterFile[MAXLEN_FILENAME]; /*!< file name of parameterfile used for starting the simulation */
extern FILE *FdInfo; /*!< file handle for info.txt log-file. */
extern FILE *FdLog; /*!< file handle for log.txt log-file. */
extern FILE *FdEnergy; /*!< file handle for energy.txt log-file. */
#ifdef SYSTEMSTATISTICS
extern FILE *FdSystem;
#endif
extern FILE *FdTimings; /*!< file handle for timings.txt log-file. */
extern FILE *FdCPU; /*!< file handle for cpu.txt log-file. */
+#ifdef SIDM_ON
+extern FILE *FdSIDMScatter; /*!< file handle for scatter.txt log-files (Different one for each core). */
+#endif
#ifdef FORCETEST
extern FILE *FdForceTest; /*!< file handle for forcetest.txt log-file. */
#endif
#ifdef SFR
extern FILE *FdSfr; /*!< file handle for sfr.txt log-file. */
#endif
#ifdef CHIMIE
extern FILE *FdChimie; /*!< file handle for chimie log-file. */
#ifdef CHIMIE_STATS
extern FILE *FdChimieStatsSNs; /*!< file handle for chimie stats-file. */
extern FILE *FdChimieStatsGas; /*!< file handle for chimie stats-file. */
long long TotNumOfForces; /*!< counts total number of force computations */
long long NumForcesSinceLastDomainDecomp; /*!< count particle updates since last domain decomposition */
/* system of units */
double G; /*!< Gravity-constant in internal units */
double UnitTime_in_s; /*!< factor to convert internal time unit to seconds/h */
double UnitMass_in_g; /*!< factor to convert internal mass unit to grams/h */
double UnitVelocity_in_cm_per_s; /*!< factor to convert intqernal velocity unit to cm/sec */
double UnitLength_in_cm; /*!< factor to convert internal length unit to cm/h */
double UnitPressure_in_cgs; /*!< factor to convert internal pressure unit to cgs units (little 'h' still around!) */
double UnitDensity_in_cgs; /*!< factor to convert internal length unit to g/cm^3*h^2 */
double UnitCoolingRate_in_cgs; /*!< factor to convert internal cooling rate to cgs units */
double UnitEnergy_in_cgs; /*!< factor to convert internal energy to cgs units */
double UnitTime_in_Megayears; /*!< factor to convert internal time to megayears/h */
double GravityConstantInternal; /*!< If set to zero in the parameterfile, the internal value of the
gravitational constant is set to the Newtonian value based on the system of
units specified. Otherwise the value provided is taken as internal gravity constant G. */
/* Cosmological parameters */
double Hubble; /*!< Hubble-constant in internal units */
double Omega0; /*!< matter density in units of the critical density (at z=0)*/
double OmegaLambda; /*!< vaccum energy density relative to crictical density (at z=0) */
double OmegaBaryon; /*!< baryon density in units of the critical density (at z=0)*/
double HubbleParam; /*!< little `h', i.e. Hubble constant in units of 100 km/s/Mpc. Only needed to get absolute physical values for cooling physics */
#ifdef HUBBLE_PARAM_WHEN_NOT_IN_UNITS
double HubbleParamNIU; /*!< little `h', i.e. Hubble constant in units of 100 km/s/Mpc. Only needed to get redshift when physical units are used */
#endif
/* Code options */
int ComovingIntegrationOn; /*!< flags that comoving integration is enabled */
int PeriodicBoundariesOn; /*!< flags that periodic boundaries are enabled */
int ResubmitOn; /*!< flags that automatic resubmission of job to queue system is enabled */
int TypeOfOpeningCriterion; /*!< determines tree cell-opening criterion: 0 for Barnes-Hut, 1 for relative criterion */
int TypeOfTimestepCriterion; /*!< gives type of timestep criterion (only 0 supported right now - unlike gadget-1.1) */
int OutputListOn; /*!< flags that output times are listed in a specified file */
/* Parameters determining output frequency */
int SnapshotFileCount; /*!< number of snapshot that is written next */
double TimeBetSnapshot; /*!< simulation time interval between snapshot files */
double TimeOfFirstSnapshot; /*!< simulation time of first snapshot files */
double CpuTimeBetRestartFile; /*!< cpu-time between regularly generated restart files */
double TimeLastRestartFile; /*!< cpu-time when last restart-file was written */
double TimeBetStatistics; /*!< simulation time interval between computations of energy statistics */
double TimeLastStatistics; /*!< simulation time when the energy statistics was computed the last time */
int NumCurrentTiStep; /*!< counts the number of system steps taken up to this point */
/* Current time of the simulation, global step, and end of simulation */
double Time; /*!< current time of the simulation */
double TimeBegin; /*!< time of initial conditions of the simulation */
double TimeStep; /*!< difference between current times of previous and current timestep */
double TimeMax; /*!< marks the point of time until the simulation is to be evolved */
/* variables for organizing discrete timeline */
double Timebase_interval; /*!< factor to convert from floating point time interval to integer timeline */
int Ti_Current; /*!< current time on integer timeline */
int Ti_nextoutput; /*!< next output time on integer timeline */
#ifdef FLEXSTEPS
int PresentMinStep; /*!< If FLEXSTEPS is used, particle timesteps are chosen as multiples of the present minimum timestep. */
int PresentMaxStep; /*!< If FLEXSTEPS is used, this is the maximum timestep in timeline units, rounded down to the next power 2 division */
#endif
#ifdef PMGRID
int PM_Ti_endstep; /*!< begin of present long-range timestep */
int PM_Ti_begstep; /*!< end of present long-range timestep */
#endif
/* Placement of PM grids */
#ifdef PMGRID
double Asmth[2]; /*!< Gives the scale of the long-range/short-range split (in mesh-cells), both for the coarse and the high-res mesh */
double Rcut[2]; /*!< Gives the maximum radius for which the short-range force is evaluated with the tree (in mesh-cells), both for the coarse and the high-res mesh */
double Corner[2][3]; /*!< lower left corner of coarse and high-res PM-mesh */
double UpperCorner[2][3]; /*!< upper right corner of coarse and high-res PM-mesh */
double Xmintot[2][3]; /*!< minimum particle coordinates both for coarse and high-res PM-mesh */
double Xmaxtot[2][3]; /*!< maximum particle coordinates both for coarse and high-res PM-mesh */
double TotalMeshSize[2]; /*!< total extension of coarse and high-res PM-mesh */
#endif
/* Variables that keep track of cumulative CPU consumption */
double TimeLimitCPU; /*!< CPU time limit as defined in parameterfile */
double CPU_TreeConstruction; /*!< time spent for constructing the gravitational tree */
double CPU_TreeWalk; /*!< actual time spent for pure tree-walks */
double CPU_Gravity; /*!< cumulative time used for gravity computation (tree-algorithm only) */
double CPU_Potential; /*!< time used for computing gravitational potentials */
double CPU_Domain; /*!< cumulative time spent for domain decomposition */
double CPU_Snapshot; /*!< time used for writing snapshot files */
double CPU_Total; /*!< cumulative time spent for domain decomposition */
double CPU_CommSum; /*!< accumulated time used for communication, and for collecting partial results, in tree-gravity */
double CPU_Imbalance; /*!< cumulative time lost accross all processors as work-load imbalance in gravitational tree */
double CPU_HydCompWalk; /*!< time used for actual SPH computations, including neighbour search */
double CPU_HydCommSumm; /*!< cumulative time used for communication in SPH, and for collecting partial results */
double CPU_HydImbalance; /*!< cumulative time lost due to work-load imbalance in SPH */
double CPU_Hydro; /*!< cumulative time spent for SPH related computations */
#ifdef SFR
double CPU_StarFormation; /*!< cumulative time spent for star formation computations */
#endif
#ifdef CHIMIE
double CPU_Chimie; /*!< cumulative time spent for chimie computations */
double CPU_ChimieDensCompWalk;
double CPU_ChimieDensCommSumm;
double CPU_ChimieDensImbalance;
double CPU_ChimieDensEnsureNgb;
double CPU_ChimieCompWalk;
double CPU_ChimieCommSumm;
double CPU_ChimieImbalance;
#endif
#ifdef MULTIPHASE
double CPU_Sticky; /*!< cumulative time spent for sticky computations */
#endif
double CPU_EnsureNgb; /*!< time needed to iterate on correct neighbour numbers */
double CPU_Predict; /*!< cumulative time to drift the system forward in time, including dynamic tree updates */
double CPU_TimeLine; /*!< time used for determining new timesteps, and for organizing the timestepping, including kicks of active particles */
double CPU_PM; /*!< time used for long-range gravitational force */
double CPU_Peano; /*!< time required to establish Peano-Hilbert order */
#ifdef DETAILED_CPU_DOMAIN
double CPU_Domain_findExtend;
double CPU_Domain_determineTopTree;
double CPU_Domain_sumCost;
double CPU_Domain_findSplit;
double CPU_Domain_shiftSplit;
double CPU_Domain_countToGo;
double CPU_Domain_exchange;
#endif
#ifdef DETAILED_CPU_GRAVITY
double CPU_Gravity_TreeWalk1;
double CPU_Gravity_TreeWalk2;
double CPU_Gravity_CommSum1;
double CPU_Gravity_CommSum2;
double CPU_Gravity_Imbalance1;
double CPU_Gravity_Imbalance2;
#endif
#ifdef COOLING
double CPU_Cooling;
#endif
#ifdef DETAILED_CPU
double CPU_Leapfrog;
double CPU_Physics;
double CPU_Residual;
double CPU_Accel;
double CPU_Begrun;
#endif
/* tree code opening criterion */
double ErrTolTheta; /*!< BH tree opening angle */
double ErrTolForceAcc; /*!< parameter for relative opening criterion in tree walk */
/* adjusts accuracy of time-integration */
double ErrTolIntAccuracy; /*!< accuracy tolerance parameter \f$ \eta \f$ for timestep criterion. The
timestep is \f$ \Delta t = \sqrt{\frac{2 \eta eps}{a}} \f$ */
double MinSizeTimestep; /*!< minimum allowed timestep. Normally, the simulation terminates if the
timestep determined by the timestep criteria falls below this limit. */
double MaxSizeTimestep; /*!< maximum allowed timestep */
double MaxRMSDisplacementFac; /*!< this determines a global timestep criterion for cosmological simulations
in comoving coordinates. To this end, the code computes the rms velocity
of all particles, and limits the timestep such that the rms displacement
is a fraction of the mean particle separation (determined from the
particle mass and the cosmological parameters). This parameter specifies
this fraction. */
double CourantFac; /*!< SPH-Courant factor */
/* frequency of tree reconstruction/domain decomposition */
double TreeDomainUpdateFrequency; /*!< controls frequency of domain decompositions */
/* Gravitational and hydrodynamical softening lengths (given in terms of an `equivalent' Plummer softening length).
* Five groups of particles are supported 0="gas", 1="halo", 2="disk", 3="bulge", 4="stars", 5="bndry"
*/
double MinGasHsmlFractional; /*!< minimum allowed SPH smoothing length in units of SPH gravitational softening length */
*GravDataIn, /*!< holds particle data to be exported to other processors */
*GravDataGet, /*!< holds particle data imported from other processors */
*GravDataResult, /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
*GravDataOut; /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */
extern struct gravdata_index
{
int Task;
int Index;
int SortIndex;
}
*GravDataIndexTable; /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */
extern struct densdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
#ifdef MULTIPHASE
int Phase;
#endif
int Index;
int Task;
#ifdef ART_CONDUCTIVITY
FLOAT EnergyIntPred;
#endif
#ifdef WITH_ID_IN_DENSITY
int ID;
#endif
}
*DensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensDataGet; /*!< holds imported particle data for SPH density computation */
*FOFGravDataIn, /*!< holds particle data to be exported to other processors */
*FOFGravDataGet, /*!< holds particle data imported from other processors */
*FOFGravDataResult, /*!< holds the partial results computed for imported particles. Note: We use GravDataResult = GravDataGet, such that the result replaces the imported data */
*FOFGravDataOut; /*!< holds partial results received from other processors. This will overwrite the GravDataIn array */
extern struct FOFgravdata_index
{
int Task;
int Index;
int SortIndex;
}
*FOFGravDataIndexTable; /*!< the particles to be exported are grouped by task-number. This table allows the results to be disentangled again and to be assigned to the correct particle */
#endif /* FOF */
#ifdef TESSEL
extern struct ghostdata_in
{
FLOAT Pos[3];
FLOAT rSearch;
int Index;
int Task;
}
*GhostDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*GhostDataGet; /*!< holds imported particle data for SPH density computation */
extern struct ghostdata_out
{
FLOAT Value;
}
*GhostDataResult, /*!< stores the locally computed SPH density results for imported particles */
*GhostDataPartialResult; /*!< imported partial SPH density results from other processors */
/* ghost particles */
//extern struct ghost_particle_data
//{
// FLOAT Pos[3]; /*!< particle position at its current time */
// FLOAT Mass; /*!< particle mass */
//}
// *gP;
extern int NumgPart;
#endif /* TESSEL */
#ifdef SYNCHRONIZE_NGB_TIMESTEP
extern struct SynchroinzeNgbTimestepdata_in
{
FLOAT Pos[3];
FLOAT Hsml;
int Ti_step;
int Ti_endstep;
int Index;
int Task;
#ifdef MULTIPHASE
int Phase;
#endif
}
*SynchroinzeNgbTimestepDataIn,
*SynchroinzeNgbTimestepDataGet;
#endif
#ifdef PY_INTERFACE
extern struct denssphdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
FLOAT Density;
FLOAT DhsmlDensityFactor;
int Index;
int Task;
FLOAT Observable;
}
*DensSphDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
*DensSphDataGet; /*!< holds imported particle data for SPH density computation */
extern struct denssphdata_out
{
FLOAT Rho;
FLOAT Div, Rot[3];
FLOAT DhsmlDensity;
FLOAT Ngb;
FLOAT GradObservable[3];
}
*DensSphDataResult, /*!< stores the locally computed SPH density results for imported particles */
*DensSphDataPartialResult; /*!< imported partial SPH density results from other processors */
extern struct sphdata_in
{
FLOAT Pos[3];
FLOAT Vel[3];
FLOAT Hsml;
FLOAT Density;
FLOAT DhsmlDensityFactor;
FLOAT ObsMoment0;
FLOAT ObsMoment1;
FLOAT Observable;
int Task;
int Index;
}
*SphDataIn, /*!< holds particle data for SPH hydro-force computation to be exported to other processors */
*SphDataGet; /*!< holds imported particle data for SPH hydro-force computation */
extern struct sphdata_out
{
FLOAT ObsMoment0;
FLOAT ObsMoment1;
FLOAT GradObservable[3];
}
*SphDataResult, /*!< stores the locally computed SPH hydro results for imported particles */
*SphDataPartialResult; /*!< imported partial SPH hydro-force results from other processors */
#endif /*PY_INTERFACE*/
#endif
diff --git a/src/begrun.c b/src/begrun.c
index eec686c..ded3882 100644
--- a/src/begrun.c
+++ b/src/begrun.c
@@ -1,2748 +1,2793 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <sys/types.h>
#include <unistd.h>
#include <gsl/gsl_rng.h>
#include "allvars.h"
#include "proto.h"
/*! \file begrun.c
* \brief initial set-up of a simulation run
*
* This file contains various functions to initialize a simulation run. In
* particular, the parameterfile is read in and parsed, the initial
* conditions or restart files are read, and global variables are
* initialized to their proper values.
*/
/*! This function performs the initial set-up of the simulation. First, the
* parameterfile is set, then routines for setting units, reading
* ICs/restart-files are called, auxialiary memory is allocated, etc.
*/
void begrun(void)
{
struct global_data_all_processes all;
#ifdef DETAILED_CPU
double tstart,tend;
tstart = second();
#endif
if(ThisTask == 0)
{
printf("\nThis is Gear, version `%s'.\n", GADGETVERSION);
printf("\nRunning on %d processors.\n", NTask);
}
read_parameter_file(ParameterFile); /* ... read in parameters for this run */
allocate_commbuffers(); /* ... allocate buffer-memory for particle
//All.StarFormationNStarsFromGas = all.StarFormationNStarsFromGas; /* do not change the param. if restarting, else, StarFormationStarMass will be wrong */
All.GrackleHSShieldingDensityThreshold = All.GrackleHSShieldingDensityThreshold*1.6726e-24/All.UnitDensity_in_cgs/(All.HubbleParam*All.HubbleParam); /* atom/cm3 to density code unit */
fprintf(FdTimings,"Step = %06d Time = %g \n\n",All.NumCurrentTiStep,All.Time);
fprintf(FdTimings,"%g %g : Total number of active particles : %d%09d\n\n",All.Time,All.TimeStep,(int) (NumActivePatricles / 1000000000), (int) (NumActivePatricles % 1000000000));