+ printf("failed to allocate memory for `StP' (%g MB).\n", bytes / (1024.0 * 1024.0));
+ 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 STELLAR_PROP
+ if(All.MaxPartStars > 0)
+ free(StP);
+#endif
+
+ if(All.MaxPartSph > 0)
+ free(SphP);
+
+ if(All.MaxPart > 0)
+ free(P);
+}
+
diff --git a/allvars.c b/allvars.c
new file mode 100644
index 0000000..895ed9e
--- /dev/null
+++ b/allvars.c
@@ -0,0 +1,329 @@
+/*! \file allvars.c
+ * \brief provides instances of all global variables.
+ */
+
+#include <stdio.h>
+#include <gsl/gsl_rng.h>
+#include "tags.h"
+#include "allvars.h"
+
+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
+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
+
+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 */
+
+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
+
+
+
+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 */
+
+
+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 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. */
+#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
+
+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 */
+
+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 */
+
+
+/* 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 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
diff --git a/allvars.h b/allvars.h
new file mode 100644
index 0000000..3d6cc5d
--- /dev/null
+++ b/allvars.h
@@ -0,0 +1,1614 @@
+
+/*! \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:
+#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 5
+#define MAXNELEMENTS 64
+#define FIRST_ELEMENT "Fe"
+#define FE 0
+#endif
+
+
+#ifdef COMPUTE_VELOCITY_DISPERSION
+#define VELOCITY_DISPERSION_SIZE 3
+#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
+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
+
+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 */
+
+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
+
+
+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 */
+}
+ *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 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. */
+#endif
+
+#ifdef MULTIPHASE
+extern FILE *FdPhase; /*!< file handle for pase.txt log-file. */
+extern FILE *FdSticky; /*!< file handle for sticky.txt log-file. */
+#endif
+
+#ifdef AGN_ACCRETION
+extern FILE *FdAccretion; /*!< file handle for accretion.txt log-file. */
+#endif
+
+#ifdef BONDI_ACCRETION
+extern FILE *FdBondi; /*!< file handle for bondi.txt log-file. */
+#endif
+
+#ifdef BUBBLES
+extern FILE *FdBubble; /*!< file handle for bubble.txt log-file. */
+#endif
+
+extern double DriftTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological drift factors */
+extern double GravKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for gravitational forces */
+extern double HydroKickTable[DRIFT_TABLE_LENGTH]; /*!< table for the cosmological kick factor for hydrodynmical forces */
+
+extern 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.
+ */
+extern struct global_data_all_processes
+{
+ long long TotNumPart; /*!< total particle numbers (global value) */
+ long long TotN_gas; /*!< total gas particle number (global value) */
+#if defined(SFR) || defined(STELLAR_PROP)
+ long long TotN_stars; /*!< total stars particle number (global value) */
+#endif
+#ifdef MULTIPHASE
+ long long TotN_sph; /*!< total sph particle number (global value) */
+ long long TotN_sticky; /*!< total sticky particle number (global value) */
+ long long TotN_stickyflaged; /*!< total sticky flaged particle number (global value) */
+ long long TotN_dark; /*!< total dark particle number (global value) */
+#endif
+
+ int MaxPart; /*!< This gives the maxmimum number of particles that can be stored on one processor. */
+ int MaxPartSph; /*!< This gives the maxmimum number of SPH particles that can be stored on one processor. */
+#ifdef STELLAR_PROP
+ int MaxPartStars; /*!< This gives the maxmimum number of Star particles that can be stored on one processor. */
+#endif
+ double BoxSize; /*!< Boxsize in case periodic boundary conditions are used */
+
+ int ICFormat; /*!< selects different versions of IC file-format */
+
+ int SnapFormat; /*!< selects different versions of snapshot file-formats */
+
+ int NumFilesPerSnapshot; /*!< number of files in multi-file snapshot dumps */
+ int NumFilesWrittenInParallel;/*!< maximum number of files that may be written simultaneously when
+ writing/reading restart-files, or when writing snapshot files */
+
+ int BufferSize; /*!< size of communication buffer in MB */
+ int BunchSizeForce; /*!< number of particles fitting into the buffer in the parallel tree-force algorithm */
+ int BunchSizeDensity; /*!< number of particles fitting into the communication buffer in the density computation */
+ int BunchSizeHydro; /*!< number of particles fitting into the communication buffer in the SPH hydrodynamical force computation */
+ int BunchSizeDomain; /*!< number of particles fitting into the communication buffer in the domain decomposition */
+#ifdef MULTIPHASE
+ int BunchSizeSticky; /*!< number of particles fitting into the communication buffer in the Chimie computation */
+#endif
+#ifdef CHIMIE
+ int BunchSizeChimie; /*!< number of particles fitting into the communication buffer in the Chimie computation */
+ int BunchSizeStarsDensity; /*!< number of particles fitting into the communication buffer in the star density computation */
+#endif
+
+
+ double PartAllocFactor; /*!< in order to maintain work-load balance, the particle load will usually
+ NOT be balanced. Each processor allocates memory for PartAllocFactor times
+ the average number of particles to allow for that */
+
+ double TreeAllocFactor; /*!< Each processor allocates a number of nodes which is TreeAllocFactor times
+ the maximum(!) number of particles. Note: A typical local tree for N
+ particles needs usually about ~0.65*N nodes. */
+
+ /* some SPH parameters */
+
+ double DesNumNgb; /*!< Desired number of SPH neighbours */
+ double MaxNumNgbDeviation; /*!< Maximum allowed deviation neighbour number */
+
+ double ArtBulkViscConst; /*!< Sets the parameter \f$\alpha\f$ of the artificial viscosity */
+ double InitGasTemp; /*!< may be used to set the temperature in the IC's */
+ double MinGasTemp; /*!< may be used to set a floor for the gas temperature */
+ double MinEgySpec; /*!< the minimum allowed temperature expressed as energy per unit mass */
+#if defined (AGN_ACCRETION) || defined (BONDI_ACCRETION)
+ double LightSpeed;
+#endif
+
+
+ /* some force counters */
+
+ 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 */
+
+
+ /* 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 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 */
+#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;
+}
+ *DensDataIn, /*!< holds particle data for SPH density computation to be exported to other processors */
+ *DensDataGet; /*!< holds imported particle data for SPH density computation */
+
+extern struct densdata_out
+{
+ FLOAT Rho;
+ FLOAT Div, Rot[3];
+ FLOAT DhsmlDensity;
+ FLOAT Ngb;
+}
+ *DensDataResult, /*!< stores the locally computed SPH density results for imported particles */
+ *DensDataPartialResult; /*!< imported partial SPH density results from other processors */
+
+
+
+extern struct hydrodata_in
+{
+ FLOAT Pos[3];
+ FLOAT Vel[3];
+ FLOAT Hsml;
+#ifdef FEEDBACK
+ FLOAT EnergySN;
+#endif
+#ifdef MULTIPHASE
+ int Phase;
+ FLOAT Entropy;
+ int StickyFlag;
+#endif
+ FLOAT Mass;
+ FLOAT Density;
+ FLOAT Pressure;
+ FLOAT F1;
+ FLOAT DhsmlDensityFactor;
+ int Timestep;
+ int Task;
+ int Index;
+#ifdef WITH_ID_IN_HYDRA
+ int ID;
+#endif
+}
+ *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 */
+
+extern struct hydrodata_out
+{
+ FLOAT Acc[3];
+ FLOAT DtEntropy;
+#ifdef FEEDBACK
+ FLOAT DtEgySpecFeedback;
+ FLOAT FeedbackAccel[3]; /*!< acceleration due to feedback force */
+ //All.StarFormationNStarsFromGas = all.StarFormationNStarsFromGas; /* do not change the param. if restarting, else, StarFormationStarMass will be wrong */
+ i = pmforce_nonperiodic(0) + pmforce_nonperiodic(1);
+ }
+ if(i != 0)
+ endrun(68688);
+#endif
+#endif
+
+
+#ifndef PERIODIC
+ if(All.ComovingIntegrationOn)
+ {
+ fac = 0.5 * All.Hubble * All.Hubble * All.Omega0;
+
+ for(i = 0; i < NumPart; i++)
+ for(j = 0; j < 3; j++)
+ P[i].GravPM[j] += fac * P[i].Pos[j];
+ }
+
+
+ /* Finally, the following factor allows a computation of cosmological simulation
+ with vacuum energy in physical coordinates */
+
+ if(All.ComovingIntegrationOn == 0)
+ {
+ fac = All.OmegaLambda * All.Hubble * All.Hubble;
+
+ for(i = 0; i < NumPart; i++)
+ for(j = 0; j < 3; j++)
+ P[i].GravPM[j] += fac * P[i].Pos[j];
+ }
+#endif
+
+}
+
+
+#endif
diff --git a/main.c b/main.c
new file mode 100644
index 0000000..2c46daa
--- /dev/null
+++ b/main.c
@@ -0,0 +1,874 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <mpi.h>
+
+#include "allvars.h"
+#include "proto.h"
+
+
+/*! \file main.c
+ * \brief start of the program
+ */
+
+/*!
+ * This function initializes the MPI communication packages, and sets
+ * cpu-time counters to 0. Then begrun() is called, which sets up
+ * the simulation either from IC's or from restart files. Finally,
+ * run() is started, the main simulation loop, which iterates over
+ * the timesteps.
+ */
+int main(int argc, char **argv)
+{
+ double t0, t1;
+
+ MPI_Init(&argc, &argv);
+ MPI_Comm_rank(MPI_COMM_WORLD, &ThisTask);
+ MPI_Comm_size(MPI_COMM_WORLD, &NTask);
+
+ if(NTask <= 1)
+ {
+ if(ThisTask == 0)
+ printf
+ ("Note: This is a massively parallel code, but you are running with 1 processor only.\nCompared to an equivalent serial code, there is some unnecessary overhead.\n");
+ }
+
+ for(PTask = 0; NTask > (1 << PTask); PTask++);
+
+ if(argc < 2)
+ {
+ if(ThisTask == 0)
+ {
+ printf("Parameters are missing.\n");
+ printf("Call with <ParameterFile> [<RestartFlag>]\n");
+ fprintf(FdLog,"%g %g : Total number of active particles : %d%09d\n\n",All.Time,All.TimeStep,(int) (NumActivePatricles / 1000000000), (int) (NumActivePatricles % 1000000000));
+#endif
+
+
+
+#ifdef DETAILED_CPU
+ tend = second();
+ All.CPU_Leapfrog += timediff(tstart, tend);
+#endif
+
+
+}
+
+
+
+/*! this function returns the next output time that is equal or larger to
+ * ti_curr
+ */
+int find_next_outputtime(int ti_curr)
+{
+
+ int i, ti, ti_next, iter = 0;
+ double next, time;
+
+ ti_next = -1;
+
+ if(All.OutputListOn)
+ {
+ for(i = 0; i < All.OutputListLength; i++)
+ {
+ time = All.OutputListTimes[i];
+
+ if(time >= All.TimeBegin && time <= All.TimeMax)
+ {
+ if(All.ComovingIntegrationOn)
+ ti = log(time / All.TimeBegin) / All.Timebase_interval;
+ else
+ ti = (time - All.TimeBegin) / All.Timebase_interval;
+
+ if(ti >= ti_curr)
+ {
+ if(ti_next == -1)
+ ti_next = ti;
+
+ if(ti_next > ti)
+ ti_next = ti;
+ }
+ }
+ }
+ }
+ else
+ {
+ if(All.ComovingIntegrationOn)
+ {
+ if(All.TimeBetSnapshot <= 1.0)
+ {
+ printf("TimeBetSnapshot > 1.0 required for your simulation.\n");
+ endrun(13123);
+ }
+ }
+ else
+ {
+ if(All.TimeBetSnapshot <= 0.0)
+ {
+ printf("TimeBetSnapshot > 0.0 required for your simulation.\n");
+ endrun(13123);
+ }
+ }
+
+ time = All.TimeOfFirstSnapshot;
+
+ iter = 0;
+
+ while(time < All.TimeBegin)
+ {
+ if(All.ComovingIntegrationOn)
+ time *= All.TimeBetSnapshot;
+ else
+ time += All.TimeBetSnapshot;
+
+ iter++;
+
+ if(iter > 1000000)
+ {
+ printf("Can't determine next output time.\n");
+ endrun(110);
+ }
+ }
+
+ while(time <= All.TimeMax)
+ {
+ if(All.ComovingIntegrationOn)
+ ti = log(time / All.TimeBegin) / All.Timebase_interval;
+ else
+ ti = (time - All.TimeBegin) / All.Timebase_interval;
+
+ if(ti >= ti_curr)
+ {
+ ti_next = ti;
+ break;
+ }
+
+ if(All.ComovingIntegrationOn)
+ time *= All.TimeBetSnapshot;
+ else
+ time += All.TimeBetSnapshot;
+
+ iter++;
+
+ if(iter > 1000000)
+ {
+ printf("Can't determine next output time.\n");
+ endrun(111);
+ }
+ }
+ }
+
+ if(ti_next == -1)
+ {
+ ti_next = 2 * TIMEBASE; /* this will prevent any further output */
+
+ if(ThisTask == 0)
+ printf("\nThere is no valid time for a further snapshot file.\n");
+ }
+ else
+ {
+ if(All.ComovingIntegrationOn)
+ next = All.TimeBegin * exp(ti_next * All.Timebase_interval);
+ else
+ next = All.TimeBegin + ti_next * All.Timebase_interval;
+
+ if(ThisTask == 0)
+ printf("\nSetting next time for snapshot file to Time_next= %g\n\n", next);
+ }
+
+ return ti_next;
+}
+
+
+
+
+/*! This routine writes one line for every timestep to two log-files. In
+ * FdInfo, we just list the timesteps that have been done, while in FdCPU the
+ * cumulative cpu-time consumption in various parts of the code is stored.
+ //printf("(%d) new star (a) i=%d id=%d (StPi=%d PIdx=%d)\n",ThisTask,N_gas+N_stars+i,P[N_gas+N_stars+i].ID,P[N_gas+N_stars+i].StPIdx,StP[P[N_gas+N_stars+i].StPIdx].PIdx);
+#endif
+
+ }
+ }
+ else /* move other particles outwards */
+ {
+
+ for(i=N_gas+N_stars;i<NumPart;i++)
+ {
+ Psave = P[i];
+ P[i] = P[i+NumNewStars];
+ P[i+NumNewStars] = Psave;
+#ifdef STELLAR_PROP
+ /* set new index for StP */
+ StP[P[i].StPIdx].PIdx = i;
+#ifdef CHECK_ID_CORRESPONDENCE
+ StP[P[i].StPIdx].ID = P[i].ID;
+#endif
+
+#endif
+ //printf("(%d) new star (b) i=%d id=%d\n",ThisTask,N_gas+N_stars+i,P[N_gas+N_stars+i].ID);