Page MenuHomec4science
Files F83534304

bondi_accretion.c
No OneTemporary
Actions

File Metadata

Created
Tue, Sep 17, 16:15

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

c4science · Help