Homec4science
Diffusion Gear (master) Browse

Gear/src/starformation.c
d524ad073023master

Last Change
Raw File

starformation.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include "allvars.h"
#include "proto.h"
#ifdef PY_INTERFACE
#include <Python.h>
#include <numpy/arrayobject.h>
#endif
#ifdef SFR
#ifdef FEEDBACK_WIND
void feedbackwind_compute_energy_kin(int mode)
{
int i;
double DeltaEgyKin;
double Tot_DeltaEgyKin;
DeltaEgyKin = 0;
Tot_DeltaEgyKin = 0;
if (mode==1)
{
LocalSysState.EnergyKin1 = 0;
LocalSysState.EnergyKin2 = 0;
}
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (mode==1)
LocalSysState.EnergyKin1 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
else
LocalSysState.EnergyKin2 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
}
}
if (mode==2)
{
DeltaEgyKin = LocalSysState.EnergyKin2 - LocalSysState.EnergyKin1;
MPI_Reduce(&DeltaEgyKin, &Tot_DeltaEgyKin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
LocalSysState.EnergyFeedbackWind -= DeltaEgyKin;
}
}
#endif
/*! \file starformation.c
* \brief Compute jet
*
*/
static double hubble_a,a3inv=1.;
/*! Compute starformation
*
*/
void star_formation(void)
{
int i,j;
int NumNewStars=0,firststar;
int *AllNumNewStars;
double T;
double p;
double dt=0;
double tstar=1.;
double MeanBaryonicDensity;
int flagST=1,flagToStar=0;
double mstar;
double mstars=0,Tot_mstars=0;
double mnewstars,Tot_mnewstars=0;
int Nnewstars,Tot_Nnewstars=0;
double star_formation_rate;
int *numlist;
double t0,t1;
double DivVel;
struct particle_data Psave;
int k;
t0 = second();
if(ThisTask == 0)
{
printf("start star formation... \n");
fflush(stdout);
}
AllNumNewStars = malloc(sizeof(int) * NTask);
#ifdef FEEDBACK_WIND
double phi,costh,sinth;
#endif
if(All.ComovingIntegrationOn)
{
hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
hubble_a = All.Hubble * sqrt(hubble_a);
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
{
hubble_a = 1.;
a3inv = 1.;
}
switch (All.StarFormationType)
{
case 0:
All.ThresholdDensity = All.StarFormationDensity;
break;
case 1:
MeanBaryonicDensity = All.OmegaBaryon * (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
All.ThresholdDensity = dmax(All.StarFormationDensity/a3inv,200*MeanBaryonicDensity);
break;
case 2:
All.ThresholdDensity = All.StarFormationDensity;
All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
break;
case 3:
All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
/* All.ThresholdDensity is fixed further */
break;
}
flagST=1;
mnewstars=0;
Nnewstars=0;
#ifdef FEEDBACK_WIND
feedbackwind_compute_energy_kin(1);
#endif
for(i = 0; i < N_gas; i++)
{
/* only active particles */
if(P[i].Ti_endstep == All.Ti_Current)
{
#ifdef CHIMIE_KINETIC_FEEDBACK
if (SphP[i].WindTime < (All.Time-All.ChimieWindTime)) /* not a wind particle */
{
#endif
/* if the particle is a false gas (new star) */
if(P[i].Type != 0)
continue;
#ifdef FEEDBACK
/* if the particle is a SN */
if(SphP[i].EnergySN)
{
/* the particle is not allowed to becomes a stars, as it is already a SN */
flagST = 0;
SphP[i].EnergySNrem = SphP[i].EnergySNrem - SphP[i].EnergySN;
/* dt in time unit here, even for ComovingIntegrationOn */
dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;
if (All.SupernovaTime==0)
SphP[i].EnergySN = SphP[i].EnergySNrem;
else
SphP[i].EnergySN = (P[i].Mass*All.SupernovaEgySpecPerMassUnit)/All.SupernovaTime*dt;
/* limit the feedback energy to the remaining energy */
if (SphP[i].EnergySN>SphP[i].EnergySNrem)
SphP[i].EnergySN = SphP[i].EnergySNrem;
if (SphP[i].EnergySNrem == 0) /* not enough energy, make a star */
{
//printf("Particle %d becomes a star\n",P[i].ID);
P[i].Type = ST;
RearrangeParticlesFlag=1;
}
}
#endif
/* check if the particle may become a star */
#ifdef SFR_NEG_DIV_ONLY
if(All.ComovingIntegrationOn)
DivVel = SphP[i].DivVel/(All.Time*All.Time) + 3*hubble_a;
else
DivVel = SphP[i].DivVel;
if (( DivVel < 0) && flagST)
#endif
{
#ifdef MULTIPHASE
if (SphP[i].Phase == GAS_STICKY)
{
#endif
if (All.StarFormationType==3)
All.ThresholdDensity = 0.5* ( pow(SphP[i].DivVel,2) + pow(SphP[i].CurlVel,2) )/(All.G);
if (SphP[i].Density*a3inv > All.ThresholdDensity)
{
/* check temperature */
#ifdef ISOTHERM_EQS
T = All.InitGasTemp;
#else
#ifdef MULTIPHASE
T = GAMMA_MINUS1*All.mumh/All.Boltzmann * SphP[i].EntropyPred;
#else
T = All.mumh/All.Boltzmann * SphP[i].EntropyPred * pow(SphP[i].Density*a3inv,GAMMA_MINUS1);
#endif
#endif
if (T < All.StarFormationTemperature)
{
/* dt in time unit here, even for ComovingIntegrationOn */
/* dt has units of Gyr/h also tstar */
dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;
tstar = All.StarFormationTime / pow(SphP[i].Density*a3inv/All.StarFormationDensity,0.5);
/* check the mass of the particle */
flagToStar = 0;
if (All.StarFormationNStarsFromGas==1)
{
/* transform the gas particle into a star particle */
mstar = P[i].Mass;
flagToStar=1;
}
else
{
if (P[i].Mass<=(1.0+All.StarFormationMgMsFraction)*All.StarFormationStarMass)
/* transform the gas particle into a star particle */
{
mstar = P[i].Mass;
flagToStar=1;
}
else
/* extract a star particle from the gas particle */
{
mstar = All.StarFormationStarMass;
flagToStar=0;
}
}
/* compute probability of forming star */
p = (P[i].Mass/mstar)*(1. - exp(-dt/tstar) );
if (get_StarFormation_random_number(P[i].ID) < p)
{
#ifdef FEEDBACK
SphP[i].EnergySNrem = P[i].Mass*All.SupernovaEgySpecPerMassUnit;
if (All.SupernovaTime==0)
SphP[i].EnergySN = SphP[i].EnergySNrem;
else
SphP[i].EnergySN = (P[i].Mass*All.SupernovaEgySpecPerMassUnit)/All.SupernovaTime*dt;
/* limit the feedback energy to the remaining energy */
if (SphP[i].EnergySN>SphP[i].EnergySNrem)
SphP[i].EnergySN = SphP[i].EnergySNrem;
if (SphP[i].EnergySNrem==0)
{
P[i].Type = ST;
RearrangeParticlesFlag=1;
}
#else
if (flagToStar) /* turns the stellar particule into star particle */
{
#ifdef STELLAR_PROP
if (N_stars+1 > All.MaxPartStars)
{
printf("No memory left for new star particle : N_stars+1=%d MaxPartStars=%d\n\n",N_stars+1,All.MaxPartStars);
endrun(999008);
}
#else
if (N_stars+1 > All.MaxPart-N_gas)
{
printf("No memory left for new star particle : N_stars+1=%d All.MaxPart-N_gas=%d\n\n",N_stars+1,All.MaxPart-N_gas);
endrun(999009);
}
#endif
P[i].Type = ST;
RearrangeParticlesFlag=1;
//printf("(%d) move gas particle to star (i=%d,id=%d)\n",ThisTask,i,P[i].ID);
/* count mass */
mnewstars+=P[i].Mass;
Nnewstars+=1;
/* gather energy */
#ifdef FEEDBACK
LocalSysState.StarEnergyFeedback += SphP[i].EgySpecFeedback * P[i].Mass;
#endif
}
else /* create a new star particle */
{
#ifdef STELLAR_PROP
if (N_stars+1 > All.MaxPartStars)
{
printf("No memory left for new star particle : N_stars+1=%d MaxPartStars=%d\n\n",N_stars+1,All.MaxPartStars);
endrun(999000);
}
#else
if (N_stars+1 > All.MaxPart-N_gas)
{
printf("No memory left for new star particle : N_stars+1=%d All.MaxPart-N_gas=%d\n\n",N_stars+1,All.MaxPart-N_gas);
endrun(999001);
}
#endif
if (NumPart+1 > All.MaxPart)
{
printf("No memory left for new particle : NumPart=%d MaxPart=%d",NumPart,All.MaxPart);
endrun(999002);
}
if (P[i].Mass <= mstar)
{
printf("Mass too small %g %g !",P[i].Mass,mstar);
endrun(999003);
}
/* the new particle get all properties of the SPH particle, except the mass */
/* the type is defined further */
j = NumPart+NumNewStars;
P[j] = P[i];
P[j].Mass = mstar;
P[j].Type = ST;
/* count mass */
mnewstars+=mstar;
Nnewstars+=1;
/* this is needed to kick parent nodes dynamically in timestep.c */
/* could be avoided if All.NumForcesSinceLastDomainDecomp > All.TotNumPart * All.TreeDomainUpdateFrequency*/
Father[j] = Father[i];
#ifdef STELLAR_PROP
/* index to Stp */
P[j].StPIdx = N_stars+NumNewStars;
/* record info at the end of StP */
/* record time */
StP[P[j].StPIdx].FormationTime = All.Time;
/* record time */
StP[P[j].StPIdx].IDProj = P[i].ID;
/* record hsml */
StP[P[j].StPIdx].Hsml = SphP[i].Hsml;
#ifdef CHIMIE
/* record initial mass */
StP[P[j].StPIdx].InitialMass = P[j].Mass;
/* record density */
StP[P[j].StPIdx].Density = SphP[i].Density;
/* record metalicity */
for(k = 0; k < NELEMENTS; k++)
StP[P[j].StPIdx].Metal[k] = SphP[i].Metal[k];
#ifdef CHIMIE_OPTIMAL_SAMPLING
init_optimal_imf_sampling(j);
#endif
#endif /* CHIMIE */
/* index to P */
StP[P[j].StPIdx].PIdx = j;
#endif
/* change proj properties */
P[i].Mass = P[i].Mass-mstar;
/* gather energy */
#ifdef FEEDBACK
LocalSysState.StarEnergyFeedback += SphP[i].EgySpecFeedback * P[j].Mass;
#endif
/* here, we should add the self force, but it is negligible */
NumNewStars++;
}
#endif
}
#ifdef FEEDBACK_WIND
else
{
tstar = tstar/All.SupernovaWindParameter;
/* here, we devide by p, because there is (1-p) prob. to be here */
p = (P[i].Mass/mstar)*(1. - exp(-dt/tstar))/(1-p);
if (get_FeedbackWind_random_number(P[i].ID) < p)
{
phi = PI*2.*get_FeedbackWind_random_number(P[i].ID+1);
costh = 1.-2.*get_FeedbackWind_random_number(P[i].ID+2);
sinth = sqrt(1.-costh*costh);
SphP[i].FeedbackWindVel[0] = All.SupernovaWindSpeed *sinth*cos(phi);
SphP[i].FeedbackWindVel[1] = All.SupernovaWindSpeed *sinth*sin(phi);
SphP[i].FeedbackWindVel[2] = All.SupernovaWindSpeed *costh;
//printf("(%d) %d vx vy vz %g %g %g %g\n",ThisTask,P[i].ID,P[i].Vel[0],P[i].Vel[1],P[i].Vel[2],sqrt( P[i].Vel[0]*P[i].Vel[0] + P[i].Vel[1]*P[i].Vel[1] + P[i].Vel[2]*P[i].Vel[2] ));
//printf("(%d) SupernovaWindSpeed %g dv=%g\n",ThisTask,All.SupernovaWindSpeed,sqrt( SphP[i].FeedbackWindVel[0]*SphP[i].FeedbackWindVel[0] + SphP[i].FeedbackWindVel[1]*SphP[i].FeedbackWindVel[1] + SphP[i].FeedbackWindVel[2]*SphP[i].FeedbackWindVel[2] ));
/* new velocities */
for(k = 0; k < 3; k++)
{
SphP[i].VelPred[k] += SphP[i].FeedbackWindVel[k];
P[i].Vel[k] += SphP[i].FeedbackWindVel[k];
}
//printf("(%d) %d vx vy vz %g %g %g %g\n",ThisTask,P[i].ID,P[i].Vel[0],P[i].Vel[1],P[i].Vel[2],sqrt( P[i].Vel[0]*P[i].Vel[0] + P[i].Vel[1]*P[i].Vel[1] + P[i].Vel[2]*P[i].Vel[2] ));
}
}
#endif
}
}
#ifdef MULTIPHASE
}
#endif
}
#ifdef CHIMIE_KINETIC_FEEDBACK
}
#endif
}
}
if (All.StarFormationNStarsFromGas != 1)
{
/* get NumNewStars from all other procs */
MPI_Allgather(&NumNewStars, 1, MPI_INT, AllNumNewStars, 1, MPI_INT, MPI_COMM_WORLD);
/* find the id of the fist local star */
firststar = All.MaxID;
for (i=0;i<ThisTask;i++)
firststar += AllNumNewStars[i];
/* now, set ids */
for (i=0;i< NumNewStars;i++)
P[NumPart+i].ID = firststar + i;
/* move the new stars to the right place
in both cases, the relati ve position of the new stars is preseved,
so, there is no need to rearrange StP.
*/
if (NumNewStars < NumPart-N_gas-N_stars) /* move stars inwards */
{
for(i=0;i<NumNewStars;i++)
{
Psave = P[N_gas+N_stars+i];
P[N_gas+N_stars+i] = P[NumPart+i];
P[NumPart+i]=Psave;
#ifdef STELLAR_PROP
/* set index for StP */
StP[P[N_gas+N_stars+i].StPIdx].PIdx = N_gas+N_stars+i;
#ifdef CHECK_ID_CORRESPONDENCE
StP[P[N_gas+N_stars+i].StPIdx].ID = P[N_gas+N_stars+i].ID;
#endif
//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 */
{
if (NumPart-N_gas-N_stars>0) /* there is other particles */
{
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);
}
/* now, set new index for the StP that where untouched */
for(i=N_gas+N_stars+1;i<N_gas+N_stars+NumNewStars;i++) /* here, we should move up to <N_gas+N_stars+NumNewStars */
{ /* however, if there is only one particle, it will not be counted */
#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
}
}
else /* there is no other particles outside, i.e, new particles are already the last ones */
{
for(i=N_gas+N_stars;i<NumPart+NumNewStars;i++)
{
#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
}
}
}
/* set new NumPart */
NumPart = NumPart+NumNewStars;
N_stars = N_stars+NumNewStars;
/*MPI_Allreduce(&NumPart, &All.TotNumPart, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);*/
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumPart, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotNumPart = 0; i < NTask; i++)
All.TotNumPart += numlist[i];
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);
}
/* compute star mass */
mstars=0;
for (i=0;i< NumPart;i++)
{
if (P[i].Type==ST)
mstars+= P[i].Mass;
}
#ifdef FEEDBACK_WIND
feedbackwind_compute_energy_kin(2);
#endif
/* share results */
MPI_Allreduce(&mstars, &Tot_mstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&mnewstars, &Tot_mnewstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&Nnewstars, &Tot_Nnewstars, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (All.StarFormationNStarsFromGas != 1)
{
#ifndef LONGIDS
All.MaxID += (unsigned int)Tot_Nnewstars;
#else
All.MaxID += (unsigned long long)Tot_Nnewstars;
#endif
}
if (All.TimeStep>0)
star_formation_rate = Tot_mnewstars/All.TimeStep;
else
star_formation_rate = 0.;
#ifndef PY_INTERFACE
if (ThisTask==0)
{
//fprintf(FdSfr, "Step %d, Time: %g, NewStarsPart: %g \n", All.NumCurrentTiStep, All.Time, Tot_mstars);
fprintf(FdSfr, "%15g %15g %15g %15g %15g\n", All.Time,All.TimeStep,Tot_mstars,Tot_mnewstars,star_formation_rate);
fflush(FdSfr);
}
#endif
//printf("(%d) N_gas = %d N_stars = %d NumPart = %d\n",ThisTask,N_gas,N_stars,NumPart);
//fflush(stdout);
if(ThisTask == 0)
{
printf("Total number of new star particles : %d \n",(int) Tot_Nnewstars);
printf("Total number of star particles : %d%09d\n",(int) (All.TotN_stars / 1000000000), (int) (All.TotN_stars % 1000000000));
fflush(stdout);
}
#ifdef CHECK_BLOCK_ORDER
int typenow;
rearrange_particle_sequence();
typenow = 0;
for(i = 0; i < NumPart; i++)
{
if ((P[i].Type<typenow)&&(typenow<=ST))
{
printf("\nBlock order error\n");
printf("(%d) i=%d id=%d Type=%d typenow=%d\n\n",ThisTask,i,P[i].ID,P[i].Type,typenow);
endrun(999004);
}
typenow = P[i].Type;
}
#endif
#ifdef CHECK_ID_CORRESPONDENCE
rearrange_particle_sequence();
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 starformation) 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(999005);
}
if(StP[P[i].StPIdx].ID != P[i].ID)
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in starformation) 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(999006);
}
}
#endif
free(AllNumNewStars);
if(ThisTask == 0)
{
printf("star formation done. \n");
fflush(stdout);
}
t1 = second();
All.CPU_StarFormation += timediff(t0, t1);
}
/*
* This routine rearrange particle sequence in order to group particles of same type
*/
void rearrange_particle_sequence(void)
{
int i,j;
struct particle_data Psave;
struct sph_particle_data Psave_sph;
int nstars=0;
int *numlist;
int sumflag;
#ifdef CHIMIE
int k;
#endif
MPI_Allreduce(&RearrangeParticlesFlag, &sumflag, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if(RearrangeParticlesFlag)
{
if(ThisTask == 0)
{
printf("Start rearrange particle sequence...\n");
fflush(stdout);
}
/* loop over all gasous particles */
for(i = 0; i < N_gas; i++)
{
if (P[i].Type == ST)
{
/* find the first non star (flag) particle */
for(j=N_gas-1;j>0;j--)
{
if (P[j].Type==0) /* we have a gasous particle */
break;
if (j==i) /* the particule is the last gaseous particle */
break; /* the particle will be inverted with itself */
if (P[j].Type==ST) /* the particule is selected to becomes a star, turn it to star */
{
#ifdef STELLAR_PROP
P[j].StPIdx = N_stars;
StP[P[j].StPIdx].FormationTime = All.Time;
StP[P[j].StPIdx].IDProj = P[j].ID;
StP[P[j].StPIdx].Hsml = SphP[j].Hsml;
#ifdef CHIMIE
/* record initial mass */
StP[P[j].StPIdx].InitialMass = P[j].Mass;
/* record density */
StP[P[j].StPIdx].Density = SphP[j].Density;
/* record metalicity */
for(k = 0; k < NELEMENTS; k++)
StP[P[j].StPIdx].Metal[k] = SphP[j].Metal[k];
#ifdef FOF_SFR
k = P[j].FoFidx;
StP[P[j].StPIdx].MassMax = FoF_P[k].MassMax;
StP[P[j].StPIdx].MassMin = FoF_P[k].MassMin;
StP[P[j].StPIdx].MassSSP = FoF_P[k].MassSSP;
StP[P[j].StPIdx].NStars = FoF_P[k].NStars;
#ifdef CHIMIE_OPTIMAL_SAMPLING
StP[P[j].StPIdx].OptIMF_k = FoF_P[k].k;
StP[P[j].StPIdx].OptIMF_CurrentMass = FoF_P[k].MassMax;
#endif
#endif
#endif /* CHIMIE */
#ifdef CHECK_ID_CORRESPONDENCE
StP[P[j].StPIdx].ID = P[j].ID;
#endif
StP[P[j].StPIdx].PIdx = j;
//printf("(%d) move to stars (b) i=%d -> j=%d id=%d Stpi=%d (N_gas=%d N_stars=%d)\n",ThisTask,j,j,P[j].ID,P[j].StPIdx,N_gas,N_stars);
#endif
nstars++;
N_gas--;
N_stars++;
continue;
}
}
//printf("(%d) --> i=%d id=%d Type=%d (N_gas=%d)\n",ThisTask,i,P[i].ID,P[i].Type,N_gas);
//printf("(%d) <-- i=%d id=%d Type=%d\n",ThisTask,j,P[j].ID,P[j].Type);
/* move the particle */
Psave = P[i]; /* copy the particle */
Psave_sph = SphP[i]; /* copy the particle */
P[i] = P[j]; /* replace by the last gaseous one */
SphP[i] = SphP[j]; /* replace by the last gaseous one */
P[j] = Psave;
#ifdef STELLAR_PROP
P[j].StPIdx = N_stars;
StP[P[j].StPIdx].FormationTime = All.Time;
StP[P[j].StPIdx].IDProj = P[j].ID;
StP[P[j].StPIdx].Hsml = Psave_sph.Hsml;
#ifdef CHIMIE
/* record initial mass */
StP[P[j].StPIdx].InitialMass = P[j].Mass;
/* record density */
StP[P[j].StPIdx].Density = Psave_sph.Density;
/* record metalicity */
for(k = 0; k < NELEMENTS; k++)
StP[P[j].StPIdx].Metal[k] = Psave_sph.Metal[k];
#ifdef FOF_SFR
k = P[j].FoFidx;
StP[P[j].StPIdx].MassMax = FoF_P[k].MassMax;
StP[P[j].StPIdx].MassMin = FoF_P[k].MassMin;
StP[P[j].StPIdx].MassSSP = FoF_P[k].MassSSP;
StP[P[j].StPIdx].NStars = FoF_P[k].NStars;
#ifdef CHIMIE_OPTIMAL_SAMPLING
StP[P[j].StPIdx].OptIMF_k= FoF_P[k].k;
StP[P[j].StPIdx].OptIMF_CurrentMass = FoF_P[k].MassMax;
#endif
#endif
#endif /* CHIMIE */
#ifdef CHECK_ID_CORRESPONDENCE
StP[P[j].StPIdx].ID = P[j].ID;
#endif
StP[P[j].StPIdx].PIdx = j;
//printf("(%d) move to stars (a) i=%d -> j=%d id=%d Stpi=%d (N_gas=%d N_stars=%d)\n",ThisTask,i,j,P[j].ID,P[j].StPIdx,N_gas,N_stars);
#endif
nstars++;
N_gas--;
N_stars++;
}
}
RearrangeParticlesFlag=0;
//if(ThisTask == 0)
// {
// printf("%d new star(s)\n",nstars);
// printf("Start rearrange particle sequence done.\n");
// fflush(stdout);
// }
}
if (sumflag) /* do this only if at least one Task has rearrange particles */
{
numlist = malloc(NTask * sizeof(int) * NTask);
/* update all.TotN_gas */
MPI_Allgather(&N_gas, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_gas = 0; i < NTask; i++)
All.TotN_gas += numlist[i];
/* update all.TotN_stars */
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);
}
}
void sfr_compute_energy_int(int mode)
{
int i;
double DeltaEgyInt;
double Tot_DeltaEgyInt;
double egyspec;
DeltaEgyInt = 0;
Tot_DeltaEgyInt = 0;
if (mode==1)
{
LocalSysState.EnergyInt1 = 0;
LocalSysState.EnergyInt2 = 0;
}
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
#ifdef MULTIPHASE
if (SphP[i].Phase== GAS_SPH)
egyspec = SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
else
egyspec = SphP[i].Entropy;
#else
egyspec = SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density * a3inv, GAMMA_MINUS1);
#endif
if (mode==1)
LocalSysState.EnergyInt1 += P[i].Mass * egyspec;
else
LocalSysState.EnergyInt2 += P[i].Mass * egyspec;
}
}
if (mode==2)
{
DeltaEgyInt = LocalSysState.EnergyInt2 - LocalSysState.EnergyInt1;
MPI_Reduce(&DeltaEgyInt, &Tot_DeltaEgyInt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
LocalSysState.StarEnergyInt -= DeltaEgyInt;
}
}
void sfr_check_number_of_stars(int mode)
{
int i;
int nstars;
int npotstars;
for(i = 0, nstars = 0, npotstars = 0; i < NumPart; i++)
{
if (P[i].Type==ST)
if (i<N_gas)
{
npotstars++;
}
else
{
nstars++;
}
}
printf("(%d) N_stars=%d npotstars+nstars=%d npotstars=%d nstars=%d\n",ThisTask,N_stars,npotstars+nstars,npotstars,nstars);
if (nstars != N_stars)
endrun(987987);
}
/****************************************************************************************/
/*
/*
/*
/* PYTHON INTERFACE
/*
/*
/*
/****************************************************************************************/
/*********************************/
/* */
/*********************************/
#ifdef PY_INTERFACE
PyObject * sfr_SetParameters(PyObject *self, PyObject *args)
{
PyObject *dict;
PyObject *key;
PyObject *value;
int ivalue;
float fvalue;
double dvalue;
/* here, we can have either arguments or dict directly */
if(PyDict_Check(args))
{
dict = args;
}
else
{
if (! PyArg_ParseTuple(args, "O",&dict))
return NULL;
}
/* check that it is a PyDictObject */
if(!PyDict_Check(dict))
{
PyErr_SetString(PyExc_AttributeError, "argument is not a dictionary.");
return NULL;
}
Py_ssize_t pos=0;
while(PyDict_Next(dict,&pos,&key,&value))
{
if(PyString_Check(key))
{
if(strcmp(PyString_AsString(key), "SfrFile")==0)
{
if(PyString_Check(value))
strcpy( All.SfrFile , PyString_AsString(value));
}
if(strcmp(PyString_AsString(key), "StarsAllocFactor")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarsAllocFactor = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationNStarsFromGas")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationNStarsFromGas = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "StarFormationMgMsFraction")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationMgMsFraction = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationStarMass")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationStarMass = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationType")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationType = PyInt_AsLong(value);
}
if(strcmp(PyString_AsString(key), "StarFormationCstar")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationCstar = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationTime")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationTime = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationDensity")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationDensity = PyFloat_AsDouble(value);
}
if(strcmp(PyString_AsString(key), "StarFormationTemperature")==0)
{
if(PyInt_Check(value)||PyLong_Check(value)||PyFloat_Check(value))
All.StarFormationTemperature = PyFloat_AsDouble(value);
}
}
}
return Py_BuildValue("i",1);
}
PyObject * sfr_GetParameters()
{
PyObject *dict;
PyObject *key;
PyObject *value;
dict = PyDict_New();
key = PyString_FromString("SfrFile");
value = PyString_FromString(All.SfrFile);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarsAllocFactor");
value = PyFloat_FromDouble(All.StarsAllocFactor);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationNStarsFromGas");
value = PyInt_FromLong(All.StarFormationNStarsFromGas);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationMgMsFraction");
value = PyFloat_FromDouble(All.StarFormationMgMsFraction);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationStarMass");
value = PyFloat_FromDouble(All.StarFormationStarMass);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationType");
value = PyInt_FromLong(All.StarFormationType);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationCstar");
value = PyFloat_FromDouble(All.StarFormationCstar);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationTime");
value = PyFloat_FromDouble(All.StarFormationTime);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationDensity");
value = PyFloat_FromDouble(All.StarFormationDensity);
PyDict_SetItem(dict,key,value);
key = PyString_FromString("StarFormationTemperature");
value = PyFloat_FromDouble(All.StarFormationTemperature);
PyDict_SetItem(dict,key,value);
return Py_BuildValue("O",dict);
}
/* definition of the method table */
static PyMethodDef sfrMethods[] = {
{"SetParameters", sfr_SetParameters, METH_VARARGS,
"set parameters"},
{"GetParameters", sfr_GetParameters, METH_VARARGS,
"get parameters"},
{NULL, NULL, 0, NULL} /* Sentinel */
};
void initsfr(void)
{
(void) Py_InitModule("sfr", sfrMethods);
import_array();
}
#endif /* PY_INTERFACE */
#endif /* SFR */
c4science · Help