Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F77925607
init.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, 06:26
Size
18 KB
Mime Type
text/x-c
Expires
Mon, Aug 19, 06:26 (2 d)
Engine
blob
Format
Raw Data
Handle
19949963
Attached To
rGEAR Gear
init.c
View Options
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file init.c
* \brief Code for initialisation of a simulation from initial conditions
*/
/*! This function reads the initial conditions, and allocates storage for the
* tree. Various variables of the particle data are initialised and An intial
* domain decomposition is performed. If SPH particles are present, the inial
* SPH smoothing lengths are determined.
*/
void
init
(
void
)
{
int
i
,
j
;
double
a3
;
#ifdef SFR
double
Mgas
=
0
,
sum_Mgas
=
0
;
int
nstars
=
0
;
int
*
numlist
;
#ifndef LONGIDS
unsigned
int
MaxID
=
0
;
#else
unsigned
long
long
MaxID
=
0
;
#endif
#endif
All
.
Time
=
All
.
TimeBegin
;
switch
(
All
.
ICFormat
)
{
case
1
:
#if (MAKEGLASS > 1)
seed_glass
();
#else
read_ic
(
All
.
InitCondFile
);
#endif
break
;
case
2
:
case
3
:
read_ic
(
All
.
InitCondFile
);
break
;
default:
if
(
ThisTask
==
0
)
printf
(
"ICFormat=%d not supported.
\n
"
,
All
.
ICFormat
);
endrun
(
0
);
}
All
.
Time
=
All
.
TimeBegin
;
All
.
Ti_Current
=
0
;
if
(
All
.
ComovingIntegrationOn
)
{
All
.
Timebase_interval
=
(
log
(
All
.
TimeMax
)
-
log
(
All
.
TimeBegin
))
/
TIMEBASE
;
a3
=
All
.
Time
*
All
.
Time
*
All
.
Time
;
}
else
{
All
.
Timebase_interval
=
(
All
.
TimeMax
-
All
.
TimeBegin
)
/
TIMEBASE
;
a3
=
1
;
}
if
(
ThisTask
==
0
)
printf
(
"
\n
Minimum Time Step (Timebase_interval) = %g
\n\n
"
,
All
.
Timebase_interval
);
set_softenings
();
All
.
NumCurrentTiStep
=
0
;
/* setup some counters */
All
.
SnapshotFileCount
=
0
;
if
(
RestartFlag
==
2
)
All
.
SnapshotFileCount
=
atoi
(
All
.
InitCondFile
+
strlen
(
All
.
InitCondFile
)
-
3
)
+
1
;
#ifdef FOF
All
.
FoF_SnapshotFileCount
=
0
;
#endif
All
.
TotNumOfForces
=
0
;
All
.
NumForcesSinceLastDomainDecomp
=
0
;
if
(
All
.
ComovingIntegrationOn
)
if
(
All
.
PeriodicBoundariesOn
==
1
)
check_omega
();
All
.
TimeLastStatistics
=
All
.
TimeBegin
-
All
.
TimeBetStatistics
;
#ifdef AGN_ACCRETION
All
.
TimeLastAccretion
=
All
.
TimeBegin
-
All
.
TimeBetAccretion
;
All
.
LastMTotInRa
=
0
;
#endif
#ifdef BONDI_ACCRETION
All
.
BondiTimeLast
=
All
.
TimeBegin
-
All
.
BondiTimeBet
;
#endif
#ifdef FOF
All
.
FoF_TimeLastFoF
=
All
.
TimeBegin
-
All
.
FoF_TimeBetFoF
;
#endif
#ifdef BUBBLES
All
.
EnergyBubbles
=
0
;
#endif
if
(
All
.
ComovingIntegrationOn
)
/* change to new velocity variable */
{
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
for
(
j
=
0
;
j
<
3
;
j
++
)
P
[
i
].
Vel
[
j
]
*=
sqrt
(
All
.
Time
)
*
All
.
Time
;
}
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
/* start-up initialization */
{
for
(
j
=
0
;
j
<
3
;
j
++
)
P
[
i
].
GravAccel
[
j
]
=
0
;
#ifdef PMGRID
for
(
j
=
0
;
j
<
3
;
j
++
)
P
[
i
].
GravPM
[
j
]
=
0
;
#endif
P
[
i
].
Ti_endstep
=
0
;
P
[
i
].
Ti_begstep
=
0
;
P
[
i
].
OldAcc
=
0
;
P
[
i
].
GravCost
=
1
;
P
[
i
].
Potential
=
0
;
#ifdef PARTICLE_FLAG
P
[
i
].
Flag
=
0
;
#endif
#ifdef TESSEL
P
[
i
].
iPref
=
-
1
;
/* index of the reference particle : -1 for normal particles */
#endif
#ifdef VANISHING_PARTICLES
P
[
i
].
VanishingFlag
=
0
;
#endif
#ifdef SFR
if
(
P
[
i
].
ID
>
MaxID
)
MaxID
=
P
[
i
].
ID
;
#endif
}
#ifdef PMGRID
All
.
PM_Ti_endstep
=
All
.
PM_Ti_begstep
=
0
;
#endif
#ifdef FLEXSTEPS
All
.
PresentMinStep
=
TIMEBASE
;
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
/* start-up initialization */
{
P
[
i
].
FlexStepGrp
=
(
int
)
(
TIMEBASE
*
get_random_number
(
P
[
i
].
ID
));
}
#endif
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
/* initialize sph_properties */
{
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
VelPred
[
j
]
=
P
[
i
].
Vel
[
j
];
SphP
[
i
].
HydroAccel
[
j
]
=
0
;
#ifdef DISSIPATION_FORCES
SphP
[
i
].
DissipationForcesAccel
[
j
]
=
0
;
#endif
}
SphP
[
i
].
DtEntropy
=
0
;
#ifdef COOLING
//SphP[i].EntropyRad = 0;
SphP
[
i
].
DtEntropyRad
=
0
;
SphP
[
i
].
DtEnergyRad
=
0
;
#endif
#ifdef DISSIPATION_FORCES
SphP
[
i
].
EnergyDissipationForces
=
0
;
SphP
[
i
].
DtEnergyDissipationForces
=
0
;
#endif
#ifdef STELLAR_FLUX
SphP
[
i
].
EnergyFlux
=
0.
;
#endif
#ifdef AGN_HEATING
SphP
[
i
].
EgySpecAGNHeat
=
0.
;
SphP
[
i
].
DtEgySpecAGNHeat
=
0.
;
#endif
#ifdef MULTIPHASE
#ifdef COUNT_COLLISIONS
SphP
[
i
].
StickyCollisionNumber
=
0
;
#endif
#endif
#ifdef FEEDBACK
SphP
[
i
].
EgySpecFeedback
=
0.
;
SphP
[
i
].
DtEgySpecFeedback
=
0.
;
SphP
[
i
].
EnergySN
=
0.
;
SphP
[
i
].
EnergySNrem
=
0.
;
SphP
[
i
].
TimeSN
=
0.
;
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
FeedbackVel
[
j
]
=
0
;
}
#endif
#ifdef FEEDBACK_WIND
for
(
j
=
0
;
j
<
3
;
j
++
)
{
SphP
[
i
].
FeedbackWindVel
[
j
]
=
0
;
}
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
SphP
[
i
].
ArtBulkViscConst
=
All
.
ArtBulkViscConstMin
;
#ifdef ART_VISCO_CD
SphP
[
i
].
ArtBulkViscConst
=
0.
;
SphP
[
i
].
DiVelAccurate
=
0.
;
SphP
[
i
].
DiVelTemp
=
0.
;
#endif
#endif
#ifdef OUTPUTOPTVAR1
SphP
[
i
].
OptVar1
=
0.
;
#endif
#ifdef OUTPUTOPTVAR2
SphP
[
i
].
OptVar2
=
0.
;
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
for
(
j
=
0
;
j
<
VELOCITY_DISPERSION_SIZE
;
j
++
)
SphP
[
i
].
VelocityDispersion
[
j
]
=
0.
;
#endif
if
(
RestartFlag
==
0
)
{
SphP
[
i
].
Hsml
=
0
;
SphP
[
i
].
Density
=
-
1
;
}
#ifdef MULTIPHASE
/* here, we set the Phase, according to the SpecificEnergy and not Entropy */
if
(
SphP
[
i
].
Entropy
>
All
.
CriticalEgySpec
)
SphP
[
i
].
Phase
=
GAS_SPH
;
/* warmer phase */
else
{
if
(
SphP
[
i
].
Entropy
>=
All
.
CriticalNonCollisionalEgySpec
)
SphP
[
i
].
Phase
=
GAS_STICKY
;
else
SphP
[
i
].
Phase
=
GAS_DARK
;
}
SphP
[
i
].
StickyFlag
=
0
;
SphP
[
i
].
StickyTime
=
All
.
Time
;
//SphP[i].StickyTime = All.Time + All.StickyIdleTime*get_random_number(P[i].ID);
#endif
#ifdef SFR
Mgas
+=
P
[
i
].
Mass
;
#endif
}
#ifdef SFR
#ifndef LONGIDS
MPI_Allreduce
(
&
MaxID
,
&
All
.
MaxID
,
1
,
MPI_INT
,
MPI_MAX
,
MPI_COMM_WORLD
);
#else
printf
(
"MPI_Allreduce must be adapted...
\n
"
);
endrun
(
-
77
);
#endif
All
.
MaxID
++
;
RearrangeParticlesFlag
=
0
;
if
(
All
.
StarFormationStarMass
==
0
)
{
/* compute the mean gas mass */
MPI_Allreduce
(
&
Mgas
,
&
sum_Mgas
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
MPI_COMM_WORLD
);
All
.
StarFormationStarMass
=
(
sum_Mgas
/
All
.
TotN_gas
)
/
All
.
StarFormationNStarsFromGas
;
}
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
/* initialize st_properties */
{
if
(
P
[
i
].
Type
==
ST
)
nstars
++
;
#ifdef STELLAR_PROP
if
(
P
[
i
].
Type
==
ST
)
{
if
(
RestartFlag
==
0
)
/* only if starting from scratch */
{
#ifndef CHIMIE_INPUT_ALL
P
[
i
].
StPIdx
=
i
-
N_gas
;
StP
[
P
[
i
].
StPIdx
].
FormationTime
=
0
;
/* bad */
StP
[
P
[
i
].
StPIdx
].
InitialMass
=
P
[
i
].
Mass
;
/* bad */
StP
[
P
[
i
].
StPIdx
].
IDProj
=
P
[
i
].
ID
;
StP
[
P
[
i
].
StPIdx
].
Hsml
=
0
;
StP
[
P
[
i
].
StPIdx
].
Density
=
-
1
;
for
(
j
=
0
;
j
<
NELEMENTS
;
j
++
)
StP
[
P
[
i
].
StPIdx
].
Metal
[
j
]
=
0.
;
StP
[
P
[
i
].
StPIdx
].
Flag
=
0
;
/*obsolete*/
#else
/* here, we restart for a file already processed by gadget */
P
[
i
].
StPIdx
=
i
-
N_gas
;
StP
[
P
[
i
].
StPIdx
].
Flag
=
0
;
/*obsolete*/
#endif
/* CHIMIE_INPUT_ALL */
}
if
(
RestartFlag
==
2
)
/* start from snapshot */
{
P
[
i
].
StPIdx
=
i
-
N_gas
;
StP
[
P
[
i
].
StPIdx
].
Flag
=
0
;
/*obsolete*/
}
StP
[
P
[
i
].
StPIdx
].
PIdx
=
i
;
#ifdef CHECK_ID_CORRESPONDENCE
StP
[
P
[
i
].
StPIdx
].
ID
=
P
[
i
].
ID
;
#endif
}
//else
// P[i].StPIdx = -1; /* shoud be set, however, may be a problem in domain.c --> must be corrected */
#endif
// STELLAR_PROP
#ifdef CHIMIE
if
(
P
[
i
].
Type
==
0
)
{
if
(
RestartFlag
==
0
&&
header
.
flag_metals
==
0
)
/* only if starting from scratch and metal block not present */
{
for
(
j
=
0
;
j
<
NELEMENTS
;
j
++
)
{
#ifdef CHIMIE_SMOOTH_METALS
SphP
[
i
].
MassMetal
[
j
]
=
(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
get_SolarMassAbundance
(
j
);
#else
SphP
[
i
].
Metal
[
j
]
=
(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
get_SolarMassAbundance
(
j
);
#endif
//if (j==FE)
// SphP[i].Metal[j] = (pow(10,All.InitGasMetallicity)-1e-10)*All.CoolingParameters_FeHSolar;
//else
// SphP[i].Metal[j] = 0;
}
}
#ifdef CHIMIE_SMOOTH_METALS
if
(
RestartFlag
==
0
&&
header
.
flag_metals
>
0
)
for
(
j
=
0
;
j
<
NELEMENTS
;
j
++
)
SphP
[
i
].
MassMetal
[
j
]
=
SphP
[
i
].
Metal
[
j
];
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
SphP
[
i
].
DeltaEgySpec
=
0
;
SphP
[
i
].
NumberOfSNIa
=
0
;
SphP
[
i
].
NumberOfSNII
=
0
;
SphP
[
i
].
SNIaThermalTime
=
-
1
;
SphP
[
i
].
SNIIThermalTime
=
-
1
;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
SphP
[
i
].
WindTime
=
All
.
TimeBegin
-
2
*
All
.
ChimieWindTime
;
SphP
[
i
].
WindFlag
=
0
;
#endif
}
#endif
/* CHIMIE */
#ifdef COOLING
#ifndef CHIMIE
All
.
GasMetal
=
(
pow
(
10
,
All
.
InitGasMetallicity
)
-
1e-10
)
*
All
.
CoolingParameters_FeHSolar
;
#endif
/* CHIMIE*/
#endif
/* COOLING */
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
if
(
P
[
i
].
Type
==
0
)
{
for
(
j
=
0
;
j
<
3
;
j
++
)
SphP
[
i
].
FeedbackUpdatedAccel
[
j
]
=
0
;
}
#endif
}
#ifdef CHECK_ID_CORRESPONDENCE
#ifdef CHIMIE
for
(
i
=
N_gas
;
i
<
N_gas
+
N_stars
;
i
++
)
{
if
(
StP
[
P
[
i
].
StPIdx
].
PIdx
!=
i
)
{
printf
(
"
\n
P/StP correspondance error
\n
"
);
printf
(
"(%d) (in domain before) 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
(
333001
);
}
if
(
StP
[
P
[
i
].
StPIdx
].
ID
!=
P
[
i
].
ID
)
{
printf
(
"
\n
P/StP correspondance error
\n
"
);
printf
(
"(%d) (in domain before) 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
(
333002
);
}
}
if
(
ThisTask
==
0
)
printf
(
"Check id correspondence before decomposition done...
\n
"
);
#endif
// CHIMIE
#endif
// CHECK_ID_CORRESPONDENCE
/* here, we would like to reduce N_stars to TotN_stars */
/* MPI_Allreduce(&N_stars, &All.TotN_stars, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); does not works */
numlist
=
malloc
(
NTask
*
sizeof
(
int
)
*
NTask
);
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
);
if
(
ThisTask
==
0
)
{
printf
(
"Total number of star particles : %d%09d
\n\n
"
,(
int
)
(
All
.
TotN_stars
/
1000000000
),
(
int
)
(
All
.
TotN_stars
%
1000000000
));
fflush
(
stdout
);
}
#endif
/*SFR*/
ngb_treeallocate
(
MAX_NGB
);
force_treeallocate
(
All
.
TreeAllocFactor
*
All
.
MaxPart
,
All
.
MaxPart
);
All
.
NumForcesSinceLastDomainDecomp
=
1
+
All
.
TotNumPart
*
All
.
TreeDomainUpdateFrequency
;
Flag_FullStep
=
1
;
/* to ensure that Peano-Hilber order is done */
domain_Decomposition
();
/* do initial domain decomposition (gives equal numbers of particles) */
ngb_treebuild
();
/* will build tree */
setup_smoothinglengths
();
#ifdef CHIMIE
#ifndef CHIMIE_INPUT_ALL
stars_setup_smoothinglengths
();
#endif
#endif
#ifdef TESSEL
setup_searching_radius
();
#endif
TreeReconstructFlag
=
1
;
/* at this point, the entropy variable normally contains the
* internal energy, read in from the initial conditions file, unless the file
* explicitly signals that the initial conditions contain the entropy directly.
* Once the density has been computed, we can convert thermal energy to entropy.
*/
#ifndef DENSITY_INDEPENDENT_SPH
/* in this case, entropy is correctely defined in setup_smoothinglengths */
#ifndef ISOTHERM_EQS
if
(
header
.
flag_entropy_instead_u
==
0
)
{
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
#ifdef MULTIPHASE
{
switch
(
SphP
[
i
].
Phase
)
{
case
GAS_SPH
:
SphP
[
i
].
Entropy
=
GAMMA_MINUS1
*
SphP
[
i
].
Entropy
/
pow
(
SphP
[
i
].
Density
/
a3
,
GAMMA_MINUS1
);
break
;
case
GAS_STICKY
:
break
;
case
GAS_DARK
:
SphP
[
i
].
Entropy
=
-
SphP
[
i
].
Entropy
;
break
;
}
}
#else
SphP
[
i
].
Entropy
=
GAMMA_MINUS1
*
SphP
[
i
].
Entropy
/
pow
(
SphP
[
i
].
Density
/
a3
,
GAMMA_MINUS1
);
#endif
}
}
#endif
#endif
/* DENSITY_INDEPENDENT_SPH */
#ifdef ENTROPYPRED
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
SphP
[
i
].
EntropyPred
=
SphP
[
i
].
Entropy
;
#endif
}
/*! This routine computes the mass content of the box and compares it to the
* specified value of Omega-matter. If discrepant, the run is terminated.
*/
void
check_omega
(
void
)
{
double
mass
=
0
,
masstot
,
omega
;
int
i
;
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
mass
+=
P
[
i
].
Mass
;
MPI_Allreduce
(
&
mass
,
&
masstot
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
MPI_COMM_WORLD
);
omega
=
masstot
/
(
All
.
BoxSize
*
All
.
BoxSize
*
All
.
BoxSize
)
/
(
3
*
All
.
Hubble
*
All
.
Hubble
/
(
8
*
M_PI
*
All
.
G
));
if
(
fabs
(
omega
-
All
.
Omega0
)
>
1.0e-3
)
{
if
(
ThisTask
==
0
)
{
printf
(
"
\n\n
I've found something odd!
\n
"
);
printf
(
"The mass content accounts only for Omega=%g,
\n
but you specified Omega=%g in the parameterfile.
\n
"
,
omega
,
All
.
Omega0
);
printf
(
"
\n
I better stop.
\n
"
);
fflush
(
stdout
);
}
endrun
(
1
);
}
}
/*! This function is used to find an initial smoothing length for each SPH
* particle. It guarantees that the number of neighbours will be between
* desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
* of the smoothing length is provided to the function density(), which will
* then iterate if needed to find the right smoothing length.
*/
void
setup_smoothinglengths
(
void
)
{
int
i
,
j
,
no
,
p
;
double
a3
;
if
(
All
.
ComovingIntegrationOn
)
{
a3
=
All
.
Time
*
All
.
Time
*
All
.
Time
;
}
else
{
a3
=
1
;
}
if
(
RestartFlag
==
0
)
{
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
{
no
=
Father
[
i
];
while
(
10
*
All
.
DesNumNgb
*
P
[
i
].
Mass
>
Nodes
[
no
].
u
.
d
.
mass
)
{
p
=
Nodes
[
no
].
u
.
d
.
father
;
if
(
p
<
0
)
break
;
no
=
p
;
}
#ifndef TWODIMS
SphP
[
i
].
Hsml
=
pow
(
3.0
/
(
4
*
M_PI
)
*
All
.
DesNumNgb
*
P
[
i
].
Mass
/
Nodes
[
no
].
u
.
d
.
mass
,
1.0
/
3
)
*
Nodes
[
no
].
len
;
#else
SphP
[
i
].
Hsml
=
pow
(
1.0
/
(
M_PI
)
*
All
.
DesNumNgb
*
P
[
i
].
Mass
/
Nodes
[
no
].
u
.
d
.
mass
,
1.0
/
2
)
*
Nodes
[
no
].
len
;
#endif
}
}
#ifdef DENSITY_INDEPENDENT_SPH
/* initialization of the entropy variable is a little trickier in this version of SPH,
since we need to make sure it 'talks to' the density appropriately */
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
SphP
[
i
].
EntVarPred
=
pow
(
SphP
[
i
].
Entropy
,
1
/
GAMMA
);
#endif
density
(
0
);
#ifdef DENSITY_INDEPENDENT_SPH
if
(
header
.
flag_entropy_instead_u
==
0
)
{
for
(
j
=
0
;
j
<
5
;
j
++
)
{
/* since ICs give energies, not entropies, need to iterate get this initialized correctly */
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
SphP
[
i
].
EntVarPred
=
pow
(
GAMMA_MINUS1
*
SphP
[
i
].
Entropy
/
pow
(
SphP
[
i
].
EgyWtDensity
/
a3
,
GAMMA_MINUS1
)
,
1
/
GAMMA
);
density
(
0
);
}
/* convert energy into entropy */
for
(
i
=
0
;
i
<
N_gas
;
i
++
)
SphP
[
i
].
Entropy
=
GAMMA_MINUS1
*
SphP
[
i
].
Entropy
/
pow
(
SphP
[
i
].
EgyWtDensity
/
a3
,
GAMMA_MINUS1
);
}
#endif
}
#ifdef CHIMIE
/*! This function is used to find an initial smoothing length for each SPH
* particle. It guarantees that the number of neighbours will be between
* desired_ngb-MAXDEV and desired_ngb+MAXDEV. For simplicity, a first guess
* of the smoothing length is provided to the function density(), which will
* then iterate if needed to find the right smoothing length.
*/
void
stars_setup_smoothinglengths
(
void
)
{
int
i
,
no
,
p
;
if
(
RestartFlag
==
0
)
{
for
(
i
=
0
;
i
<
NumPart
;
i
++
)
{
if
(
P
[
i
].
Type
==
ST
)
{
no
=
Father
[
i
];
while
(
10
*
All
.
DesNumNgb
*
P
[
i
].
Mass
>
Nodes
[
no
].
u
.
d
.
mass
)
{
p
=
Nodes
[
no
].
u
.
d
.
father
;
if
(
p
<
0
)
break
;
no
=
p
;
}
#ifndef TWODIMS
StP
[
P
[
i
].
StPIdx
].
Hsml
=
pow
(
3.0
/
(
4
*
M_PI
)
*
All
.
DesNumNgb
*
P
[
i
].
Mass
/
Nodes
[
no
].
u
.
d
.
mass
,
1.0
/
3
)
*
Nodes
[
no
].
len
;
#else
StP
[
P
[
i
].
StPIdx
].
Hsml
=
pow
(
1.0
/
(
M_PI
)
*
All
.
DesNumNgb
*
P
[
i
].
Mass
/
Nodes
[
no
].
u
.
d
.
mass
,
1.0
/
2
)
*
Nodes
[
no
].
len
;
#endif
}
}
}
stars_density
();
}
#endif
/*! If the code is run in glass-making mode, this function populates the
* simulation box with a Poisson sample of particles.
*/
#if (MAKEGLASS > 1)
void
seed_glass
(
void
)
{
int
i
,
k
,
n_for_this_task
;
double
Range
[
3
],
LowerBound
[
3
];
double
drandom
,
partmass
;
long
long
IDstart
;
All
.
TotNumPart
=
MAKEGLASS
;
partmass
=
All
.
Omega0
*
(
3
*
All
.
Hubble
*
All
.
Hubble
/
(
8
*
M_PI
*
All
.
G
))
*
(
All
.
BoxSize
*
All
.
BoxSize
*
All
.
BoxSize
)
/
All
.
TotNumPart
;
All
.
MaxPart
=
All
.
PartAllocFactor
*
(
All
.
TotNumPart
/
NTask
);
/* sets the maximum number of particles that may */
allocate_memory
();
header
.
npartTotal
[
1
]
=
All
.
TotNumPart
;
header
.
mass
[
1
]
=
partmass
;
if
(
ThisTask
==
0
)
{
printf
(
"
\n
Glass initialising
\n
PartMass= %g
\n
"
,
partmass
);
printf
(
"TotNumPart= %d%09d
\n\n
"
,
(
int
)
(
All
.
TotNumPart
/
1000000000
),
(
int
)
(
All
.
TotNumPart
%
1000000000
));
}
/* set the number of particles assigned locally to this task */
n_for_this_task
=
All
.
TotNumPart
/
NTask
;
if
(
ThisTask
==
NTask
-
1
)
n_for_this_task
=
All
.
TotNumPart
-
(
NTask
-
1
)
*
n_for_this_task
;
NumPart
=
0
;
IDstart
=
1
+
(
All
.
TotNumPart
/
NTask
)
*
ThisTask
;
/* split the temporal domain into Ntask slabs in z-direction */
Range
[
0
]
=
Range
[
1
]
=
All
.
BoxSize
;
Range
[
2
]
=
All
.
BoxSize
/
NTask
;
LowerBound
[
0
]
=
LowerBound
[
1
]
=
0
;
LowerBound
[
2
]
=
ThisTask
*
Range
[
2
];
srand48
(
ThisTask
);
for
(
i
=
0
;
i
<
n_for_this_task
;
i
++
)
{
for
(
k
=
0
;
k
<
3
;
k
++
)
{
drandom
=
drand48
();
P
[
i
].
Pos
[
k
]
=
LowerBound
[
k
]
+
Range
[
k
]
*
drandom
;
P
[
i
].
Vel
[
k
]
=
0
;
}
P
[
i
].
Mass
=
partmass
;
P
[
i
].
Type
=
1
;
P
[
i
].
ID
=
IDstart
+
i
;
NumPart
++
;
}
}
#endif
Event Timeline
Log In to Comment