Homec4science
Diffusion Gear (master) Browse

Gear/src/new.phase.c
fb822919c90cmaster

Last Change
Raw File

new.phase.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 MULTIPHASE
/*! \file phase.c
* \brief Compute phase mixing
*
*/
/*! compute phase mixing
*
*/
void update_phase()
{
int i;
double a3,hubble_a;
#ifdef PHASE_MIXING
FLOAT EgySpec=0;
#endif
double dt;
#ifdef COLDGAS_CYCLE
double Pwc,Pcw;
double flux_in_cgs,X,tau;
#endif
double Pcol,ColTime;
double DensityNorm;
int *numlist;
N_sph = 0;
N_sticky = 0;
N_stickyflaged = 0;
N_dark = 0;
NumColPotLocal = 0;
if(All.ComovingIntegrationOn)
{
a3 = All.Time * All.Time * All.Time;
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);
}
else
{
a3 = 1;
hubble_a = 1;
}
for(i = 0; i < N_gas; i++)
{
//if(P[i].Ti_endstep == All.Ti_Current) /* if the particle is active */
{
SphP[i].StickyFlag = 0;
#ifdef PHASE_MIXING
#ifdef FEEDBACK
if((P[i].Type == 0) && (SphP[i].EnergySN==0))
#else
if(P[i].Type == 0)
#endif
{
switch (SphP[i].Phase)
{
case GAS_SPH:
/* compute temperature from entropy */
EgySpec = 1./GAMMA_MINUS1 * pow(SphP[i].Density / a3,GAMMA_MINUS1)* SphP[i].Entropy;
if (EgySpec < All.CriticalEgySpec) /* to sticky phase */
{
SphP[i].Phase = GAS_STICKY;
SphP[i].StickyFlag = 0;
SphP[i].StickyTime = All.Time; /* may be elligible for collision */
SphP[i].Entropy = EgySpec;
}
break;
case GAS_STICKY:
/* compute temperature from specific energy */
EgySpec = SphP[i].Entropy;
if (EgySpec >= All.CriticalEgySpec) /* to sph phase */
{
SphP[i].Phase = GAS_SPH;
SphP[i].Entropy = GAMMA_MINUS1 * EgySpec / pow(SphP[i].Density / a3, GAMMA_MINUS1);
}
break;
}
#ifdef COLDGAS_CYCLE
/*
Here, we do nothing during the first step, because the UV flux is unknown
*/
if (((SphP[i].Phase==GAS_STICKY)||(SphP[i].Phase==GAS_DARK))&&All.NumCurrentTiStep>0)
{
/* cold/warm gas cycle */
/* compute dt */
dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;
X = 0.;
#ifdef STELLAR_FLUX
/* compute stellar flux */
flux_in_cgs = SphP[i].EnergyFlux* All.UnitEnergy_in_cgs/All.UnitTime_in_s/pow(All.UnitLength_in_cm, 2);
X = X + flux_in_cgs/C / All.HeatingPeSolarEnergyDensity;
#endif
#ifdef EXTERNAL_FLUX
/* compute outer flux */
X = X + All.HeatingExternalFLuxEnergyDensity/All.HeatingPeSolarEnergyDensity ;
#endif
if (X>0)
{
/* warm gas */
tau = pow(X,+All.ColdGasCycleTransitionParameter) * All.ColdGasCycleTransitionTime;
Pwc = 1.0-exp(-dt/tau);
/* cold gas */
tau = pow(X,-All.ColdGasCycleTransitionParameter) * All.ColdGasCycleTransitionTime;
Pcw = 1.0-exp(-dt/tau);
}
else
{
Pwc = 1.0;
Pcw = 0.0;
}
/* compute transition */
switch(SphP[i].Phase)
{
case GAS_STICKY:
if (get_random_number(P[i].ID) < Pwc)
{
SphP[i].Phase = GAS_DARK;
}
break;
case GAS_DARK:
if (get_random_number(P[i].ID) < Pcw)
{
SphP[i].Phase = GAS_STICKY;
SphP[i].StickyFlag = 0;
SphP[i].StickyTime = All.Time; /* may be elligible for collision */
}
break;
}
}
#endif /* COLDGAS_CYCLE */
#endif /* PHASE_MIXING */
/* determine if a sticky particle may collide or not */
if(SphP[i].Phase==GAS_STICKY)
{
if(SphP[i].StickyTime <= All.Time) /* may collide */
{
/* new implementation */
SphP[i].StickyFlag = 1;
/* this was the old implementation */
//dt = (All.Ti_Current - P[i].Ti_begstep) * All.Timebase_interval / hubble_a;
//
///* because mean free path depends on density */
//if (All.StickyDensity>0)
// {
// //ColTime = All.StickyTime * pow(All.StickyDensity/SphP[i].Density,All.StickyDensityPower);
// DensityNorm = SphP[i].Density/All.StickyDensity;
// ColTime = All.StickyTime / (DensityNorm * pow( (1+pow(DensityNorm/3,2)) ,1/2.));
// if (DensityNorm > 10)
// {
// ColTime = ColTime*1e10; /* cut off */
// }
// }
//else
// ColTime = All.StickyTime;
//
//
//Pcol = 1.0-exp(-dt/ColTime);
//if (get_random_number(P[i].ID) < Pcol)
// {
// SphP[i].StickyFlag = 1; /* the particle may collide */
// SphP[i].StickyTime = All.Time + All.StickyIdleTime;
// NumColPotLocal++;
// }
}
}
} /* end of Ptype==0 */
} /* end of active particles */
/* count number of each type */
if(P[i].Type == 0)
{
switch (SphP[i].Phase)
{
case GAS_SPH:
N_sph++;
break;
case GAS_STICKY:
N_sticky++;
if (SphP[i].StickyFlag == 1)
N_stickyflaged++;
break;
case GAS_DARK:
N_dark++;
break;
}
}
} /* end of main loop */
/* share results */
//MPI_Allreduce(&N_sph, &All.TotN_sph, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
//MPI_Allreduce(&N_sticky, &All.TotN_sticky, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
//MPI_Allreduce(&N_stickyflaged, &All.TotN_stickyflaged, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
//MPI_Allreduce(&N_dark, &All.TotN_dark, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&N_sph, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_sph = 0; i < NTask; i++)
All.TotN_sph += numlist[i];
MPI_Allgather(&N_sticky, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_sticky = 0; i < NTask; i++)
All.TotN_sticky += numlist[i];
MPI_Allgather(&N_stickyflaged, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_stickyflaged = 0; i < NTask; i++)
All.TotN_stickyflaged += numlist[i];
MPI_Allgather(&N_dark, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, All.TotN_dark = 0; i < NTask; i++)
All.TotN_dark += numlist[i];
free(numlist);
if (ThisTask==0)
{
fprintf(FdPhase, "Step %d, Time: %g GasPart: %d SphPart: %d StickyPart: %d DarkPart: %d StickyPartFlaged: %d\n", All.NumCurrentTiStep, All.Time, (int)All.TotN_gas, (int)All.TotN_sph, (int)All.TotN_sticky,(int)All.TotN_dark, (int)All.TotN_stickyflaged);
fflush(FdPhase);
}
}
#endif
c4science · Help