Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77955043
new.phase.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Aug 17, 11:15
Size
8 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 11:15 (2 d)
Engine
blob
Format
Raw Data
Handle
19952401
Attached To
rGEAR Gear
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
Event Timeline
Log In to Comment