Page MenuHomec4science
Files F65233263

init.c
No OneTemporary
Actions

File Metadata

Created
Sun, Jun 2, 02:08

init.c
View Options

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file init.c
* \brief Code for initialisation of a simulation from initial conditions
*/
/*! This function reads the initial conditions, and allocates storage for the
* tree. Various variables of the particle data are initialised and An intial
* domain decomposition is performed. If SPH particles are present, the inial
* SPH smoothing lengths are determined.
*/
void init(void)
{
int i, j;
double a3;
#ifdef SFR
double Mgas=0,sum_Mgas=0;
int nstars=0;
int *numlist;
#ifndef LONGIDS
unsigned int MaxID=0;
#else
unsigned long long MaxID=0;
#endif
#endif
All.Time = All.TimeBegin;
switch (All.ICFormat)
{
case 1:
#if (MAKEGLASS > 1)
seed_glass();
#else
read_ic(All.InitCondFile);
#endif
break;
case 2:
case 3:
read_ic(All.InitCondFile);
break;
default:
if(ThisTask == 0)
printf("ICFormat=%d not supported.\n", All.ICFormat);
endrun(0);
}
All.Time = All.TimeBegin;
All.Ti_Current = 0;
if(All.ComovingIntegrationOn)
{
All.Timebase_interval = (log(All.TimeMax) - log(All.TimeBegin)) / TIMEBASE;
a3 = All.Time * All.Time * All.Time;
}
else
{
All.Timebase_interval = (All.TimeMax - All.TimeBegin) / TIMEBASE;
a3 = 1;
}
if (ThisTask==0)
printf("\nMinimum Time Step (Timebase_interval) = %g \n\n", All.Timebase_interval);
set_softenings();
All.NumCurrentTiStep = 0; /* setup some counters */
All.SnapshotFileCount = 0;
if(RestartFlag == 2)
All.SnapshotFileCount = atoi(All.InitCondFile + strlen(All.InitCondFile) - 3) + 1;
#ifdef FOF
All.FoF_SnapshotFileCount = 0;
#endif
All.TotNumOfForces = 0;
All.NumForcesSinceLastDomainDecomp = 0;
if(All.ComovingIntegrationOn)
if(All.PeriodicBoundariesOn == 1)
check_omega();
All.TimeLastStatistics = All.TimeBegin - All.TimeBetStatistics;
#ifdef AGN_ACCRETION
All.TimeLastAccretion = All.TimeBegin - All.TimeBetAccretion;
All.LastMTotInRa = 0;
#endif
#ifdef BONDI_ACCRETION
All.BondiTimeLast = All.TimeBegin - All.BondiTimeBet;
#endif
#ifdef CHIMIE
#ifdef CHIMIE_TIMEBET
All.ChimieTimeLastChimie = All.TimeBegin - All.ChimieTimeBetChimie; /* shoud be change following FoF_TimeLastFoF */
All.ChimieTi_begstep = All.Ti_Current;
#endif
#endif
#ifdef FOF
All.FoF_TimeLastFoF = All.TimeBegin;
#endif
#ifdef BUBBLES
All.EnergyBubbles=0;
#endif
if(All.ComovingIntegrationOn) /* change to new velocity variable */
{
for(i = 0; i < NumPart; i++)
for(j = 0; j < 3; j++)
P[i].Vel[j] *= sqrt(All.Time) * All.Time;
}
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
for(j = 0; j < 3; j++)
P[i].GravAccel[j] = 0;
#ifdef PMGRID
for(j = 0; j < 3; j++)
P[i].GravPM[j] = 0;
#endif
P[i].Ti_endstep = 0;
P[i].Ti_begstep = 0;
P[i].OldAcc = 0;
P[i].GravCost = 1;
P[i].Potential = 0;
#ifdef PARTICLE_FLAG
P[i].Flag = 0;
#endif
#ifdef TESSEL
P[i].iPref = -1; /* index of the reference particle : -1 for normal particles */
#endif
#ifdef VANISHING_PARTICLES
P[i].VanishingFlag=0;
#endif
#ifdef SFR
if (P[i].ID>MaxID)
MaxID=P[i].ID;
#endif
}
#ifdef PMGRID
All.PM_Ti_endstep = All.PM_Ti_begstep = 0;
#endif
#ifdef FLEXSTEPS
All.PresentMinStep = TIMEBASE;
for(i = 0; i < NumPart; i++) /* start-up initialization */
{
P[i].FlexStepGrp = (int) (TIMEBASE * get_random_number(P[i].ID));
}
#endif
for(i = 0; i < N_gas; i++) /* initialize sph_properties */
{
for(j = 0; j < 3; j++)
{
SphP[i].VelPred[j] = P[i].Vel[j];
SphP[i].HydroAccel[j] = 0;
#ifdef DISSIPATION_FORCES
SphP[i].DissipationForcesAccel[j] = 0;
#endif
}
SphP[i].DtEntropy = 0;
#ifdef COOLING
//SphP[i].EntropyRad = 0;
SphP[i].DtEntropyRad = 0;
SphP[i].DtEnergyRad = 0;
#endif
#ifdef DISSIPATION_FORCES
SphP[i].EnergyDissipationForces = 0;
SphP[i].DtEnergyDissipationForces = 0;
#endif
#ifdef STELLAR_FLUX
SphP[i].EnergyFlux = 0.;
#endif
#ifdef AGN_HEATING
SphP[i].EgySpecAGNHeat = 0.;
SphP[i].DtEgySpecAGNHeat = 0.;
#endif
#ifdef MULTIPHASE
#ifdef COUNT_COLLISIONS
SphP[i].StickyCollisionNumber = 0;
#endif
#endif
#ifdef FEEDBACK
SphP[i].EgySpecFeedback = 0.;
SphP[i].DtEgySpecFeedback = 0.;
SphP[i].EnergySN = 0.;
SphP[i].EnergySNrem = 0.;
SphP[i].TimeSN = 0.;
for(j = 0; j < 3; j++)
{
SphP[i].FeedbackVel[j] = 0;
}
#endif
#ifdef FEEDBACK_WIND
for(j = 0; j < 3; j++)
{
SphP[i].FeedbackWindVel[j] = 0;
}
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
SphP[i].ArtBulkViscConst = All.ArtBulkViscConstMin;
#ifdef ART_VISCO_CD
SphP[i].ArtBulkViscConst = 0.;
SphP[i].DiVelAccurate = 0.;
SphP[i].DiVelTemp = 0.;
#endif
#endif
#ifdef OUTPUTOPTVAR1
SphP[i].OptVar1 = 0.;
#endif
#ifdef OUTPUTOPTVAR2
SphP[i].OptVar2 = 0.;
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
for(j = 0; j < VELOCITY_DISPERSION_SIZE; j++)
SphP[i].VelocityDispersion[j] = 0.;
#endif
if(RestartFlag == 0)
{
SphP[i].Hsml = 0;
SphP[i].Density = -1;
}
#ifdef MULTIPHASE
/* here, we set the Phase, according to the SpecificEnergy and not Entropy */
if (SphP[i].Entropy > All.CriticalEgySpec)
SphP[i].Phase = GAS_SPH; /* warmer phase */
else
{
if (SphP[i].Entropy >= All.CriticalNonCollisionalEgySpec)
SphP[i].Phase = GAS_STICKY;
else
SphP[i].Phase = GAS_DARK;
}
SphP[i].StickyFlag = 0;
SphP[i].StickyTime = All.Time;
//SphP[i].StickyTime = All.Time + All.StickyIdleTime*get_random_number(P[i].ID);
#endif
#ifdef SFR
Mgas += P[i].Mass;
#endif
}
#ifdef SFR
#ifndef LONGIDS
MPI_Allreduce(&MaxID, &All.MaxID, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD);
#else
printf("MPI_Allreduce must be adapted...\n");
endrun(-77);
#endif
All.MaxID++;
RearrangeParticlesFlag=0;
if (All.StarFormationStarMass==0)
{
/* compute the mean gas mass */
MPI_Allreduce(&Mgas, &sum_Mgas, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
All.StarFormationStarMass = (sum_Mgas/All.TotN_gas) / All.StarFormationNStarsFromGas;
}
for(i = 0; i < NumPart; i++) /* initialize st_properties */
{
if (P[i].Type==ST)
nstars++;
#ifdef STELLAR_PROP
if (P[i].Type==ST)
{
if (RestartFlag==0) /* only if starting from scratch */
{
#ifndef CHIMIE_INPUT_ALL
P[i].StPIdx = i-N_gas;
StP[P[i].StPIdx].FormationTime = 0; /* bad */
StP[P[i].StPIdx].InitialMass = P[i].Mass; /* bad */
StP[P[i].StPIdx].IDProj = P[i].ID;
StP[P[i].StPIdx].Hsml = 0;
StP[P[i].StPIdx].Density = -1;
for(j = 0; j < NELEMENTS; j++)
StP[P[i].StPIdx].Metal[j] = 0.;
StP[P[i].StPIdx].Flag = 0; /*obsolete*/
#else /* here, we restart for a file already processed by gadget */
P[i].StPIdx = i-N_gas;
StP[P[i].StPIdx].Flag = 0; /*obsolete*/
#endif /* CHIMIE_INPUT_ALL */
}
if (RestartFlag==2) /* start from snapshot */
{
P[i].StPIdx = i-N_gas;
StP[P[i].StPIdx].Flag = 0; /*obsolete*/
}
StP[P[i].StPIdx].PIdx = i;
#ifdef CHECK_ID_CORRESPONDENCE
StP[P[i].StPIdx].ID = P[i].ID;
#endif
}
//else
// P[i].StPIdx = -1; /* shoud be set, however, may be a problem in domain.c --> must be corrected */
#endif // STELLAR_PROP
#ifdef CHIMIE
if (P[i].Type==0)
{
if (RestartFlag==0 && header.flag_metals==0) /* only if starting from scratch and metal block not present */
{
for(j = 0; j < NELEMENTS; j++)
{
#ifdef CHIMIE_SMOOTH_METALS
SphP[i].MassMetal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*get_SolarMassAbundance(j);
#else
SphP[i].Metal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*get_SolarMassAbundance(j);
#endif
//if (j==FE)
// SphP[i].Metal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar;
//else
// SphP[i].Metal[j] = 0;
}
}
if (RestartFlag==2)
{
if(SphP[i].Metal[0] == 0)
{
for(j = 0; j < NELEMENTS; j++)
{
#ifdef CHIMIE_SMOOTH_METALS
SphP[i].MassMetal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*get_SolarMassAbundance(j);
#else
SphP[i].Metal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*get_SolarMassAbundance(j);
#endif
}
}
else
{
#ifdef CHIMIE_SMOOTH_METALS
for(j = 0; j < NELEMENTS; j++)
{
SphP[i].MassMetal[j] = SphP[i].Metal[j];
}
#endif
}
}
#ifdef CHIMIE_SMOOTH_METALS
if (RestartFlag==0 && header.flag_metals>0)
for(j = 0; j < NELEMENTS; j++)
SphP[i].MassMetal[j] = SphP[i].Metal[j];
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
SphP[i].DeltaEgySpec = 0;
SphP[i].NumberOfSNIa = 0;
SphP[i].NumberOfSNII = 0;
SphP[i].SNIaThermalTime = -1;
SphP[i].SNIIThermalTime = -1;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
SphP[i].WindTime = All.TimeBegin-2*All.ChimieWindTime;
SphP[i].WindFlag = 0;
#endif
}
#endif /* CHIMIE */
#ifdef COOLING
#ifndef CHIMIE
All.GasMetal = (pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar;
#endif /* CHIMIE*/
#endif /* COOLING */
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
if (P[i].Type==0)
{
for (j=0;j<3;j++)
SphP[i].FeedbackUpdatedAccel[j] = 0;
}
#endif
}
#ifdef CHECK_ID_CORRESPONDENCE
#ifdef CHIMIE
for(i = N_gas; i < N_gas+N_stars; i++)
{
if( StP[P[i].StPIdx].PIdx != i )
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in domain before) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
endrun(333001);
}
if(StP[P[i].StPIdx].ID != P[i].ID)
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in domain before) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
endrun(333002);
}
}
if (ThisTask==0)
printf("Check id correspondence before decomposition done...\n");
#endif // CHIMIE
#endif // CHECK_ID_CORRESPONDENCE
/* here, we would like to reduce N_stars to TotN_stars */
/* MPI_Allreduce(&N_stars, &All.TotN_stars, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); does not works */
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&N_stars, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_stars = 0; i < NTask; i++)
All.TotN_stars += numlist[i];
free(numlist);
if(ThisTask == 0)
{
printf("Total number of star particles : %d%09d\n\n",(int) (All.TotN_stars / 1000000000), (int) (All.TotN_stars % 1000000000));
fflush(stdout);
}
#endif /*SFR*/
ngb_treeallocate(MAX_NGB);
force_treeallocate(All.TreeAllocFactor * All.MaxPart, All.MaxPart);
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
Flag_FullStep = 1; /* to ensure that Peano-Hilber order is done */
domain_Decomposition(); /* do initial domain decomposition (gives equal numbers of particles) */
ngb_treebuild(); /* will build tree */
setup_smoothinglengths();
#ifdef CHIMIE
#ifndef CHIMIE_INPUT_ALL
stars_setup_smoothinglengths();
#endif
#endif
#ifdef TESSEL
setup_searching_radius();
#endif
TreeReconstructFlag = 1;
/* at this point, the entropy variable normally contains the
* internal energy, read in from the initial conditions file, unless the file
* explicitly signals that the initial conditions contain the entropy directly.
* Once the density has been computed, we can convert thermal energy to entropy.
*/
#ifndef DENSITY_INDEPENDENT_SPH /* in this case, entropy is correctely defined in setup_smoothinglengths */
#ifndef ISOTHERM_EQS
if(header.flag_entropy_instead_u == 0)
{
for(i = 0; i < N_gas; i++)
{
#ifdef MULTIPHASE
{
switch(SphP[i].Phase)
{
case GAS_SPH:
SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
break;
case GAS_STICKY:
break;
case GAS_DARK:
SphP[i].Entropy = -SphP[i].Entropy;
break;
}
}
#else
SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
}
}
#endif
#endif /* DENSITY_INDEPENDENT_SPH */
#ifdef ENTROPYPRED
for(i = 0; i < N_gas; i++)
SphP[i].EntropyPred = SphP[i].Entropy;
#endif
#ifdef TRACE_ACC
trace_pos();
#endif
}
/*! This routine computes the mass content of the box and compares it to the
* specified value of Omega-matter. If discrepant, the run is terminated.
*/
void check_omega(void)
{
double mass = 0, masstot, omega;
int i;
for(i = 0; i < NumPart; i++)
mass += P[i].Mass;
MPI_Allreduce(&mass, &masstot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
omega =
masstot / (All.BoxSize * All.BoxSize * All.BoxSize) / (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
if(fabs(omega - All.Omega0) > 1.0e-3)
{
if(ThisTask == 0)
{
printf("\n\nI've found something odd!\n");
printf
("The mass content accounts only for Omega=%g,\nbut you specified Omega=%g in the parameterfile.\n",
omega, All.Omega0);
printf("\nI better stop.\n");
fflush(stdout);
}
endrun(1);
}
}
/*! This function is used to find an initial smoothing length for each SPH
* particle. It guarantees that the number of neighbours will be between
* desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
* of the smoothing length is provided to the function density(), which will
* then iterate if needed to find the right smoothing length.
*/
void setup_smoothinglengths(void)
{
int i, j, no, p;
double a3;
if(All.ComovingIntegrationOn)
{
a3 = All.Time * All.Time * All.Time;
}
else
{
a3 = 1;
}
if(RestartFlag == 0 || RestartFlag == 2)
{
for(i = 0; i < N_gas; i++)
{
no = Father[i];
while(10 * All.DesNumNgb * P[i].Mass > Nodes[no].u.d.mass)
{
p = Nodes[no].u.d.father;
if(p < 0)
break;
no = p;
}
#ifndef TWODIMS
SphP[i].Hsml =
pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
#else
SphP[i].Hsml =
pow(1.0 / (M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 2) * Nodes[no].len;
#endif
}
}
#ifdef DENSITY_INDEPENDENT_SPH
/* initialization of the entropy variable is a little trickier in this version of SPH,
since we need to make sure it 'talks to' the density appropriately */
for(i = 0; i < N_gas; i++)
SphP[i].EntVarPred = pow(SphP[i].Entropy,1/GAMMA);
#endif
density(0);
#ifdef DENSITY_INDEPENDENT_SPH
if(header.flag_entropy_instead_u == 0)
{
for(j=0;j<5;j++)
{ /* since ICs give energies, not entropies, need to iterate get this initialized correctly */
for(i = 0; i < N_gas; i++)
SphP[i].EntVarPred = pow( GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].EgyWtDensity/a3 , GAMMA_MINUS1) ,1/GAMMA);
density(0);
}
/* convert energy into entropy */
for(i = 0; i < N_gas; i++)
SphP[i].Entropy = GAMMA_MINUS1 * SphP[i].Entropy / pow(SphP[i].EgyWtDensity/a3 , GAMMA_MINUS1);
}
#endif
}
#ifdef CHIMIE
/*! This function is used to find an initial smoothing length for each SPH
* particle. It guarantees that the number of neighbours will be between
* desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
* of the smoothing length is provided to the function density(), which will
* then iterate if needed to find the right smoothing length.
*/
void stars_setup_smoothinglengths(void)
{
int i, no, p;
if(RestartFlag == 0 || RestartFlag == 2)
{
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == ST)
{
no = Father[i];
while(10 * All.DesNumNgb * P[i].Mass > Nodes[no].u.d.mass)
{
p = Nodes[no].u.d.father;
if(p < 0)
break;
no = p;
}
#ifndef TWODIMS
StP[P[i].StPIdx].Hsml =
pow(3.0 / (4 * M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 3) * Nodes[no].len;
#else
StP[P[i].StPIdx].Hsml =
pow(1.0 / (M_PI) * All.DesNumNgb * P[i].Mass / Nodes[no].u.d.mass, 1.0 / 2) * Nodes[no].len;
#endif
}
}
}
stars_density();
}
#endif
/*! If the code is run in glass-making mode, this function populates the
* simulation box with a Poisson sample of particles.
*/
#if (MAKEGLASS > 1)
void seed_glass(void)
{
int i, k, n_for_this_task;
double Range[3], LowerBound[3];
double drandom, partmass;
long long IDstart;
All.TotNumPart = MAKEGLASS;
partmass = All.Omega0 * (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G))
* (All.BoxSize * All.BoxSize * All.BoxSize) / All.TotNumPart;
All.MaxPart = All.PartAllocFactor * (All.TotNumPart / NTask); /* sets the maximum number of particles that may */
allocate_memory();
header.npartTotal[1] = All.TotNumPart;
header.mass[1] = partmass;
if(ThisTask == 0)
{
printf("\nGlass initialising\nPartMass= %g\n", partmass);
printf("TotNumPart= %d%09d\n\n",
(int) (All.TotNumPart / 1000000000), (int) (All.TotNumPart % 1000000000));
}
/* set the number of particles assigned locally to this task */
n_for_this_task = All.TotNumPart / NTask;
if(ThisTask == NTask - 1)
n_for_this_task = All.TotNumPart - (NTask - 1) * n_for_this_task;
NumPart = 0;
IDstart = 1 + (All.TotNumPart / NTask) * ThisTask;
/* split the temporal domain into Ntask slabs in z-direction */
Range[0] = Range[1] = All.BoxSize;
Range[2] = All.BoxSize / NTask;
LowerBound[0] = LowerBound[1] = 0;
LowerBound[2] = ThisTask * Range[2];
srand48(ThisTask);
for(i = 0; i < n_for_this_task; i++)
{
for(k = 0; k < 3; k++)
{
drandom = drand48();
P[i].Pos[k] = LowerBound[k] + Range[k] * drandom;
P[i].Vel[k] = 0;
}
P[i].Mass = partmass;
P[i].Type = 1;
P[i].ID = IDstart + i;
NumPart++;
}
}
#endif

Event Timeline

c4science · Help