Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65941773
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
Fri, Jun 7, 05:55
Size
8 KB
Mime Type
text/x-c
Expires
Sun, Jun 9, 05:55 (2 d)
Engine
blob
Format
Raw Data
Handle
18149016
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