Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F83525116
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
Tue, Sep 17, 15:04
Size
8 KB
Mime Type
text/x-c
Expires
Thu, Sep 19, 15:04 (2 d)
Engine
blob
Format
Raw Data
Handle
20852212
Attached To
R195 SCITAS application benchmarks
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;
int Tis,Tic;
#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 */
if(P[i].Type == 0)
{
SphP[i].StickyFlag = 0;
#ifdef PHASE_MIXING
#ifdef FEEDBACK
if(SphP[i].EnergySN==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;
/* disable :Mon Feb 21 16:47:28 CET 2011 */ //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;
/* disable :Mon Feb 21 16:47:28 CET 2011 */ //SphP[i].StickyTime = All.Time; /* may be elligible for collision */
}
break;
}
}
#endif
#endif
/* determine if a sticky particle may collide or not */
if(SphP[i].Phase==GAS_STICKY)
{
if (All.ComovingIntegrationOn)
{
Tis = log(SphP[i].StickyTime/All.TimeBegin) / All.Timebase_interval;
Tic = All.Ti_Current;
dt = get_cosmictime_difference(Tis,Tic);
}
else
{
dt = All.Time - SphP[i].StickyTime;
}
if(dt >= All.StickyIdleTime)
SphP[i].StickyFlag = 1;
// /* disable :Mon Feb 21 16:47:28 CET 2011 */ //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 active particles */ /* P[i].Type==0 */
/* 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