Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85979268
bondi_accretion.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
Thu, Oct 3, 10:11
Size
4 KB
Mime Type
text/x-c
Expires
Sat, Oct 5, 10:11 (2 d)
Engine
blob
Format
Raw Data
Handle
21279498
Attached To
rPNBODY pNbody
bondi_accretion.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 BONDI_ACCRETION
/*! \file bondi_accretion.c
* \brief Compute bondi accretion on central black hole
*
*/
static
double
hubble_a
,
a3inv
;
void
compute_bondi_accretion
()
{
int
i
,
n
,
j
;
double
pressure
,
density
;
//double pressure_max,density_max;
//double *density_list,*pressure_list,*radius_list;
double
pressure_sum
,
density_sum
;
double
BondiMdot
;
int
numngb
,
startnode
,
phase
,
ngb
,
ngb_sum
;
FLOAT
pos
[
3
],
h_i
,
h_i2
;
double
dx
,
dy
,
dz
;
double
h_j
,
r
,
r2
;
double
u
,
wk
,
hinv
,
hinv3
;
if
(
All
.
ComovingIntegrationOn
)
{
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
);
a3inv
=
1
/
(
All
.
Time
*
All
.
Time
*
All
.
Time
);
}
else
a3inv
=
hubble_a
=
1
;
if
((
All
.
Time
-
All
.
BondiTimeLast
)
>=
All
.
BondiTimeBet
||
All
.
NumCurrentTiStep
==
0
)
{
All
.
BondiTimeLast
+=
All
.
BondiTimeBet
;
density
=
0
;
pressure
=
0
;
ngb
=
0
;
/* find neighbors around the center */
pos
[
0
]
=
0
;
pos
[
1
]
=
0
;
pos
[
2
]
=
0
;
h_i
=
All
.
SofteningTable
[
0
]
*
All
.
BondiHsmlFactor
;
h_i2
=
h_i
*
h_i
;
hinv
=
1.0
/
h_i
;
hinv3
=
hinv
*
hinv
*
hinv
;
startnode
=
All
.
MaxPart
;
#ifdef MULTIPHASE
phase
=
GAS_SPH
;
numngb
=
ngb_treefind_phase_pairs
(
&
pos
[
0
],
h_i
,
phase
,
&
startnode
);
#else
phase
=
0
;
numngb
=
ngb_treefind_pairs
(
&
pos
[
0
],
h_i
,
phase
,
&
startnode
);
#endif
for
(
n
=
0
;
n
<
numngb
;
n
++
)
{
j
=
Ngblist
[
n
];
dx
=
pos
[
0
]
-
P
[
j
].
Pos
[
0
];
dy
=
pos
[
1
]
-
P
[
j
].
Pos
[
1
];
dz
=
pos
[
2
]
-
P
[
j
].
Pos
[
2
];
r2
=
dx
*
dx
+
dy
*
dy
+
dz
*
dz
;
h_j
=
SphP
[
j
].
Hsml
;
if
(
r2
<
h_i2
)
{
r
=
sqrt
(
r2
);
u
=
r
*
hinv
;
ngb
++
;
if
(
u
<
0.5
)
{
wk
=
hinv3
*
(
KERNEL_COEFF_1
+
KERNEL_COEFF_2
*
(
u
-
1
)
*
u
*
u
);
}
else
{
wk
=
hinv3
*
KERNEL_COEFF_5
*
(
1.0
-
u
)
*
(
1.0
-
u
)
*
(
1.0
-
u
);
}
density
+=
P
[
j
].
Mass
*
wk
;
pressure
+=
SphP
[
j
].
Pressure
*
P
[
j
].
Mass
/
SphP
[
j
].
Density
*
wk
;
}
}
/* sum contribution from other nodes */
MPI_Allreduce
(
&
density
,
&
density_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
pressure
,
&
pressure_sum
,
1
,
MPI_DOUBLE
,
MPI_SUM
,
MPI_COMM_WORLD
);
MPI_Allreduce
(
&
ngb
,
&
ngb_sum
,
1
,
MPI_INT
,
MPI_SUM
,
MPI_COMM_WORLD
);
density
=
density_sum
;
pressure
=
pressure_sum
;
/************************************/
/* loop over all particles */
/************************************/
/*
density_max = 0;
pressure_max = 0;
for(i = 0; i < NumPart; i++)
{
if(P[i].Type == 0)
{
if (density_max < SphP[i].Density)
{
density_max = SphP[i].Density;
pressure_max = SphP[i].Pressure;
}
}
}
density_list = malloc(NTask * sizeof(double));
pressure_list = malloc(NTask * sizeof(double));
radius_list = malloc(NTask * sizeof(double));
MPI_Allgather(&density_max, 1, MPI_DOUBLE, density_list, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Allgather(&pressure_max, 1, MPI_DOUBLE, pressure_list, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Allgather(&radius_max, 1, MPI_DOUBLE, radius_list, 1, MPI_DOUBLE, MPI_COMM_WORLD);
density = 0;
pressure = 0;
for(i = 0; i < NTask; i++)
if (density < density_list[i])
{
density = density_list[i];
pressure= pressure_list[i];
radius= radius_list[i];
}
free(density_list);
free(pressure_list);
free(radius_list);
*/
/* compute bondi accretion */
BondiMdot
=
4
*
PI
*
pow
(
All
.
G
*
All
.
BondiBlackHoleMass
,
2
)
*
density
*
pow
(
GAMMA
*
pressure
/
density
,
-
3.
/
2.
);
/* compute bondi power */
All
.
BondiPower
=
All
.
BondiEfficiency
*
(
All
.
LightSpeed
*
All
.
LightSpeed
)
*
BondiMdot
;
/* output info */
if
(
ThisTask
==
0
)
{
fprintf
(
FdBondi
,
"Step: %d, Time: %g, BondiMdot: %g, BondiPower: %g, BondiDensity: %g, BondiPressure: %g NumNgb: %d
\n
"
,
All
.
NumCurrentTiStep
,
All
.
Time
,
BondiMdot
,
All
.
BondiPower
,
density
,
pressure
,
ngb_sum
);
fflush
(
FdBondi
);
}
}
}
#endif
Event Timeline
Log In to Comment