Page MenuHomec4science
Files F77544389

agn_feedback.c
No OneTemporary
Actions

File Metadata

Created
Thu, Aug 15, 02:29

agn_feedback.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 AGN_ACCRETION
/*! \file bubble.c
* \brief Init bubble
*
*/
/*! compute bubble
*
*/
void compute_agn_accretion()
{
int i,k;
double AccretionRadius2;
double Ma[7],Ma_sum[7];
double r2,r;
#ifdef AGN_USE_ANGULAR_MOMENTUM
double Lvec[3],Lvec_sum[3];
double theta,phi;
#endif
/*double a3inv=1.;*/ /* !!!!!!!!!!!! */
if((All.Time - All.TimeLastAccretion) >= All.TimeBetAccretion || All.NumCurrentTiStep==0)
{
All.TimeLastAccretion += All.TimeBetAccretion;
if(ThisTask == 0)
printf("computing accretion\n");
AccretionRadius2 = All.AccretionRadius*All.AccretionRadius;
for (k = 0; k < 7; k++)
{
Ma[k]=0.;
Ma_sum[k]=0.;
}
#ifdef AGN_USE_ANGULAR_MOMENTUM
Lvec[0] = Lvec[1] = Lvec[2] = 0;
#endif
/*
loop over all particles
*/
for(i = 0; i < NumPart; i++)
{
r2 = P[i].Pos[0]*P[i].Pos[0] + P[i].Pos[1]*P[i].Pos[1] + P[i].Pos[2]*P[i].Pos[2];
if (r2<AccretionRadius2)
{
#ifdef MULTIPHASE
if (P[i].Type==0)
{
switch (SphP[i].Phase)
{
case GAS_SPH:
Ma[0] += P[i].Mass;
break;
case GAS_STICKY:
case GAS_DARK:
Ma[6] += P[i].Mass;
break;
}
}
else
{
Ma[P[i].Type] += P[i].Mass;
}
#else
Ma[P[i].Type] += P[i].Mass;
#endif
#ifdef AGN_USE_ANGULAR_MOMENTUM
Lvec[0] += P[i].Pos[1]*P[i].Vel[2] - P[i].Pos[2]*P[i].Vel[1];
Lvec[1] += P[i].Pos[2]*P[i].Vel[0] - P[i].Pos[0]*P[i].Vel[2];
Lvec[2] += P[i].Pos[0]*P[i].Vel[1] - P[i].Pos[1]*P[i].Vel[0];
#endif
}
}
/* share results */
MPI_Allreduce(Ma, Ma_sum, 7, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#ifdef AGN_USE_ANGULAR_MOMENTUM
MPI_Allreduce(&Lvec, &Lvec_sum, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
#endif
/*
compute parameters for feedback
All.LastMTotInRa : total Mass in Ra at last feedback
All.MTotInRa : total Mass in Ra
All.dMTotInRa : diff between total Mass in Ra now and total Mass in Ra last feedback
All.C : speed of light
All.AGNFactor : factor efficiency
All.MinMTotInRa : minimum accumulated mass in Ra needed to generate a buble
*/
/* store results Ma_sum[6]=sticky */
All.MTotInRa = Ma_sum[0]+Ma_sum[6]+Ma_sum[1]+Ma_sum[2]+Ma_sum[3]+Ma_sum[4];
if (All.NumCurrentTiStep==0)
{
All.LastMTotInRa = All.MTotInRa;
All.dMTotInRa = 0;
}
else
{
All.dMTotInRa = All.MTotInRa - All.LastMTotInRa;
}
/*
if(ThisTask == 0)
printf("--> time=%g dMTotInRa=%g MTotInRa=%g LastMTotInRa=%g MinMTotInRa=%g AGNFactor=%g\n",All.Time,All.dMTotInRa,All.MTotInRa,All.LastMTotInRa,All.MinMTotInRa,All.AGNFactor);
*/
#ifdef AGN_FEEDBACK
/* check if it is time for a bubble */
if (All.dMTotInRa > All.MinMTotInRa)
{
/* compute jet energy */
All.BubblesE[All.BubblesIndex] = All.AGNFactor * All.dMTotInRa * All.LightSpeed* All.LightSpeed;
All.LastMTotInRa = All.MTotInRa;
#ifdef AGN_USE_ANGULAR_MOMENTUM
/* compute jet direction */
r = sqrt( Lvec_sum[0]*Lvec_sum[0] + Lvec_sum[1]*Lvec_sum[1] + Lvec_sum[2]*Lvec_sum[2] );
theta = acos(Lvec_sum[2]/r);
r = sqrt( Lvec_sum[0]*Lvec_sum[0] + Lvec_sum[1]*Lvec_sum[1] );
phi = acos(Lvec_sum[0]/r);
All.BubblesA[All.BubblesIndex] = theta;
All.BubblesB[All.BubblesIndex] = phi;
#endif
/* set time for the bubbles */
All.BubblesTime[All.BubblesIndex]= 0.;
}
#endif
/* write output */
if(ThisTask == 0)
{
#ifdef MULTIPHASE
fprintf(FdAccretion, "Step: %d, Time: %g, M1a: %g, M1b: %g, M2: %g, M3: %g, M4: %g, M5: %g, M6: %g\n", All.NumCurrentTiStep, All.Time,
Ma_sum[0],Ma_sum[6],Ma_sum[1],Ma_sum[2],Ma_sum[3],Ma_sum[4],Ma_sum[5]);
#else
fprintf(FdAccretion, "Step: %d, Time: %g, M1: %g, M2: %g, M3: %g, M4: %g, M5: %g, M6: %g\n", All.NumCurrentTiStep, All.Time,
Ma_sum[0],Ma_sum[1],Ma_sum[2],Ma_sum[3],Ma_sum[4],Ma_sum[5]);
#endif
fflush(FdAccretion);
}
}
}
#endif

Event Timeline

c4science · Help