From 05bf72a9e7195d675caf902cb63e14bc57a94161 Mon Sep 17 00:00:00 2001 From: Benjamin Moster Date: Wed, 8 Apr 2020 21:27:03 +0200 Subject: [PATCH] Compute the stellar mass left in a population (MassLeft) more accurately by integrating over each timestep --- src/allvars.c | 1 + src/allvars.h | 4 ++- src/galaxies.c | 3 +++ src/read_trees.c | 65 ++++++++++++++++++++++++++++++++++++++++-------- src/setup.c | 1 + 5 files changed, 62 insertions(+), 12 deletions(-) diff --git a/src/allvars.c b/src/allvars.c index f37eb12..4e46888 100644 --- a/src/allvars.c +++ b/src/allvars.c @@ -51,6 +51,7 @@ float *CosmicTime; // Cosmic time corresponding to the above scale fact float *Timestep; // Timesteps corresponding the above scale factors float *DynTime; // Dynamical halo times corresponding the above scale factors +float *PopAge; // Age of a stellar population float *MassLeft; // Stores the fraction of the mass left between each scale factor float *MassFormed; // Stores the mass of all populations ever formed in a branch float *ICMFormed; // Stores the mass of all populations in the ICM diff --git a/src/allvars.h b/src/allvars.h index 0c1bf3b..ae5bb7a 100644 --- a/src/allvars.h +++ b/src/allvars.h @@ -36,13 +36,14 @@ #define CODE "EMERGE" ///< Code name #define BRANCH "master" ///< Branch name -#define VERSION "1.0.1" ///< Code version +#define VERSION "1.0.2" ///< Code version #define SYMBOL '#' ///< Symbol used to start output lines #define NSTRING 200 ///< Number of elements used for a generic string #define SSTRING 50 ///< Number of elements used for a short string #define LINESIZE 100 ///< Number of characters used for the screen output #define ALLOC_TOLERANCE 0.2 ///< Tolerance in the memory allocation +#define NINTQUAD 16 ///< Number of integration nodes for gaussian quadrature #define ZMAX_SMFERR 4.0 ///< Maximum redshift up to which the observed stellar mass error increases #define SSFRTHRESH 0.3 ///< Fraction of the inverse Hubble time used for the SSFR threshold @@ -542,6 +543,7 @@ extern float *CosmicTime; ///< Cosmic time corresponding to the above s extern float *Timestep; ///< Timesteps corresponding the above scale factors extern float *DynTime; ///< Dynamical halo times corresponding the above scale factors +extern float *PopAge; ///< Age of a stellar population extern float *MassLeft; ///< Stores the fraction of the mass left between each scale factor extern float *MassFormed; ///< Stores the mass of all populations ever formed in a branch extern float *ICMFormed; ///< Stores the mass of all populations in the ICM diff --git a/src/galaxies.c b/src/galaxies.c index aeaf0fc..ec9cb2b 100644 --- a/src/galaxies.c +++ b/src/galaxies.c @@ -46,6 +46,9 @@ void init_galaxies(void) int iforest; #endif + //Print what is done... + if (ThisTask == 0) printf("%s\n%s Populating dark matter haloes with galaxies...\n", All.fullline,All.startline); + //First find the maximum number of leaves #ifndef COMPUTE_ICM //Loop through all trees diff --git a/src/read_trees.c b/src/read_trees.c index 3a549b9..ca3c045 100644 --- a/src/read_trees.c +++ b/src/read_trees.c @@ -21,6 +21,7 @@ #include #include #include +#include #include #include "allvars.h" #include "proto.h" @@ -731,7 +732,7 @@ void empty_read_buffer(int offset, int hc) H[offset + n].mmp = h[n].mmp; H[offset + n].a = h[n].a; - H[offset + n].mvir = h[n].mvir/All.h_100/All.m_unit;; + H[offset + n].mvir = h[n].mvir/All.h_100/All.m_unit; H[offset + n].rvir = h[n].rvir/All.h_100/All.x_unit/1000.0*h[n].a; H[offset + n].c = h[n].c; H[offset + n].lambda = h[n].lambda; @@ -1471,7 +1472,11 @@ void sort_by_size(void) //And increment the number of forests one last time iforest++; - if (iforest != Nforests) endrun("The number of forests after sorting is not the same as the stored number"); + if (iforest != Nforests) + { + printf("%d %d %d\n",ThisTask,iforest,Nforests); + endrun("The number of forests after sorting is not the same as the stored number"); + } #endif //Finally restore all the values that have been used for the sorting (icoprog, idescgal and possibly imass) @@ -1897,9 +1902,13 @@ void distribute_trees(int ntask, double load) void get_timesteps(void) { - int i, j, nTimesteps, maxnTimesteps; + int i, j, m, nTimesteps, maxnTimesteps; + float t0, t1, ml; float *a, *alist, *alistall; + double *xxx, *www; char buf[500]; + gsl_integration_fixed_workspace *w; + const gsl_integration_fixed_type *T = gsl_integration_fixed_legendre; //Print what is done... if (ThisTask == 0) printf("%s\n%s Identifying the number of timesteps and the list of scale factors from the trees...\n", All.fullline,All.startline); @@ -2066,22 +2075,57 @@ void get_timesteps(void) } #endif + //Compute the population age table which stores the age of each stellar population + PopAge = emalloc_movable(&PopAge, "PopAge", All.NTimesteps * All.NTimesteps * sizeof(float)); + for (i = 0; i < All.NTimesteps; i++) + { //The age at the first snapshot is the age of the Universe + PopAge[i*All.NTimesteps] = log10(CosmicTime[i]*All.t_unit); + for (j = 1; j < All.NTimesteps; j++) + { //For each other snapshot, the population age is the time since the previous snapshot + if (j <= i) PopAge[i*All.NTimesteps+j] = log10((CosmicTime[i]-CosmicTime[j-1])*All.t_unit); + else PopAge[i*All.NTimesteps+j] = 0.0; + } + } + //Compute the MassLeft table which stores the fraction of mass left after any time interval MassLeft = emalloc_movable(&MassLeft,"MassLeft", All.NTimesteps * All.NTimesteps * sizeof(float)); + //Go through all timesteps for (i = 0; i < All.NTimesteps; i++) - { + { //Go through all populations for this timestep for (j = 0; j < All.NTimesteps; j++) - { - if (j < i) MassLeft[i*All.NTimesteps+j] = 1.0 - 0.05 * log( ((CosmicTime[i]-CosmicTime[j])*All.t_unit/1.4e6) + 1.0); - else MassLeft[i*All.NTimesteps+j] = 1.0; - } - } + { //If the population forms after or at this snapshot set t0 to 0, otherwise to the next snapshot (lower bound for age) + if (i <= j) t0 = 0.0; + else t0 = pow(10.0,PopAge[i*All.NTimesteps+j+1])/All.t_unit; + //Set t1 to the current population age (upper bound for age) + t1 = pow(10.0,PopAge[i*All.NTimesteps+j])/All.t_unit; + //If population forms before or at this snapshot + if (i >= j) + { //Allocate the integration workspace and link the quadrature nodes and weights + w = gsl_integration_fixed_alloc(T, NINTQUAD, t0, t1, 0.0, 0.0); + xxx = gsl_integration_fixed_nodes(w); + www = gsl_integration_fixed_weights(w); + //Set initial fraction of mass left to 0 + ml = 0.0; + //Integrate over the quadrature nodes by multiplying the weights and f_left and summing it up + for (m = 0; m < NINTQUAD; m++) ml += www[m] * (1.0 - 0.05 * log( (xxx[m]*All.t_unit/1.4e6) + 1.0)); + //Divide by the timestep width and write the result into the MassLeft array + MassLeft[i*All.NTimesteps+j] = ml / (t1 - t0); + //Free the integration workspace + gsl_integration_fixed_free(w); + } + else + { //Otherwise set the fraction of stars left to 1 (will never be used though) + MassLeft[i*All.NTimesteps+j] = 1.0; + } + }//End j + }//End i //Broadcast so that every universe has the joy MPI_Bcast(ScaleFactor, All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(CosmicTime, All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(Timestep, All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(DynTime, All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); + MPI_Bcast(PopAge, All.NTimesteps * All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(MassLeft, All.NTimesteps * All.NTimesteps, MPI_FLOAT, 0, MPI_COMM_WORLD); MPI_Bcast(Output_iscale, All.Noutputredshifts, MPI_INT, 0, MPI_COMM_WORLD); MPI_Bcast(OutputRedshifts, All.Noutputredshifts, MPI_FLOAT, 0, MPI_COMM_WORLD); @@ -2224,7 +2268,6 @@ void add_orphans(void) else H[iorphan].descid = maxid+2; //We now use the upid of the orphan to record the true descendants halo ID H[iorphan].upid = H[idesc].haloid; - //Copy all relevant information from the true descendant H[iorphan].iscale = H[idesc].iscale; H[iorphan].type = 2; //All orphans are type 2 @@ -2526,7 +2569,7 @@ void compute_halo_history(void) //Go through all haloes for (ihalo = 0; ihalo < Nhalos; ihalo++) { - //If this halo does not have any progenitors set its iprog to -1 and its impeak to this halo's index + //If this halo does not have any progenitors set its impeak to this halo's index if (H[ihalo].np < 1) H[ihalo].impeak = ihalo; //Copy the index of the peak halo mass from the progenitor else H[ihalo].impeak = H[H[ihalo].iprog].impeak; diff --git a/src/setup.c b/src/setup.c index 114e584..8914961 100644 --- a/src/setup.c +++ b/src/setup.c @@ -803,6 +803,7 @@ void setup(void){ void finalize (void) { efree_movable(MassLeft); + efree_movable(PopAge); efree_movable(DynTime); efree_movable(Timestep); efree_movable(CosmicTime);