https://github.com/N-BodyShop/changa
Tip revision: 69a5c80002f52a61a7002dc8fb77cde41c61b97e authored by Tim Haines on 23 October 2018, 23:21:25 UTC
Complete overhaul of build system
Complete overhaul of build system
Tip revision: 69a5c80
sinks.cpp
/*
* Sink particle module originally written for GASOLINE by James
* Wadsley and Jillian Bellovary (Black holes).
* Modified for inclusion in ChaNGa
*/
#include <math.h>
#ifdef HAVE_VALUES_H
#include <values.h>
#else
#include <float.h>
#endif
#include "ParallelGravity.h"
#include "smooth.h"
#include "sinks.h"
#include "Sph.h"
#include "formatted_string.h"
/// Utility function
static inline double sqr(double x)
{
return x*x;
}
///
/// @brief initialize parameters for sinks
///
void Sinks::AddParams(PRM prm, Parameters ¶ms)
{
//
// Sink parameters
//
bDoSinks = 0;
prmAddParam(prm,"bDoSinks",paramBool,&bDoSinks,sizeof(int),
"sinks","enable/disable sinks = -sinks");
bBHSink = 0;
prmAddParam(prm,"bBHSink",paramBool,&bBHSink,sizeof(int),
"bhsink","Bondi-Hoyle type sink = -bhsink");
bSmallBHSmooth = 0;
prmAddParam(prm,"bSmallBHSmooth",paramBool,&bSmallBHSmooth,
sizeof(int), "smallbhsmooth", "smooth BH feedback over blastwave or smoothing radius = -smallbhsmooth");
bBHTurnOffCooling = 0;
prmAddParam(prm,"bBHTurnOffCooling",paramBool,&bBHTurnOffCooling,
sizeof(int), "bhturnoffcooling",
"turn off cooling for BHs = -bhturnoffcooling");
bDoBHKick = 0;
prmAddParam(prm,"bDoBHKick",paramBool,&bDoBHKick,sizeof(int),
"dobhkick","turn on grav recoil for mergers = -dobhkick");
dBHSinkEddEff = 0.1;
prmAddParam(prm,"dBHSinkEddEff",paramDouble,&dBHSinkEddEff,
sizeof(double),"bhsinkeddeff",
"<BHSink Eddington Efficiency>");
dBHSinkFeedbackEff = 0.05;
prmAddParam(prm,"dBHSinkFeedbackEff",paramDouble,
&dBHSinkFeedbackEff,sizeof(double),"bhsinkfbeff",
"<BHSink Feedback Efficiency>");
dBHSinkAlpha = 1;
prmAddParam(prm,"dBHSinkAlpha",paramDouble,&dBHSinkAlpha,
sizeof(double),"bhsinkalpha", "<BHSink Alpha>");
bBHMindv = 0;
prmAddParam(prm,"bBHMindv",paramBool,&bBHMindv,sizeof(int),
"bhmindv","use mindeltav for BH accretion = -bhmindv");
bBHAccreteAll = 1;
prmAddParam(prm,"bBHAccreteAll",paramBool,&bBHAccreteAll,
sizeof(int),"bhaccreteall",
"BHs can accrete any particle = -bhaccreteall");
bDoSinksAtStart = 0;
prmAddParam(prm,"bDoSinksAtStart",paramBool,&bDoSinksAtStart,
sizeof(int),
"sinksas","enable/disable sinks at start = -sinksas");
bSinkThermal = 0;
prmAddParam(prm,"bSinkThermal",paramBool,&bSinkThermal,
sizeof(int), "tsinks",
"enable/disable thermal energy in sink calcs = -tsinks");
dSinkRadius = 0.0;
prmAddParam(prm,"dSinkRadius",paramDouble,&dSinkRadius,
sizeof(double),"sinkr", "<Sink Radius>");
dSinkBoundOrbitRadius = 0.0;
prmAddParam(prm,"dSinkBoundOrbitRadius",paramDouble,
&dSinkBoundOrbitRadius,sizeof(double),"sinkbor",
"<Sink Bound Orbit Radius>");
dSinkMustAccreteRadius = 0.0;
prmAddParam(prm,"dSinkMustAccreteRadius",paramDouble,
&dSinkMustAccreteRadius,sizeof(double),"sinkr",
"<Sink Radius>");
dDeltaSink = params.dDelta;
prmAddParam(prm,"dDeltaSink", paramDouble, &dDeltaSink,
sizeof(double), "dDeltaSink",
"<Maximum sink timestep in years> = dDelta");
dSinkMassMin = FLT_MAX; /* Needs to be explicitly set */
prmAddParam(prm,"dSinkMassMin", paramDouble, &dSinkMassMin,
sizeof(double), "dSinkMassMin",
"<Minimum Mass to act as a sink> = FLT_MAX" );
iSinkRung = 0;
prmAddParam(prm,"iSinkRung", paramInt, &iSinkRung, sizeof(int),
"iSinkRung", "<Sink Rung> = 0");
bSinkForm = 0;
prmAddParam(prm,"bSinkForm",paramBool,&bSinkForm,sizeof(int),
"sinkform","enable/disable sinks = -sinkform");
nJeans = 50;
prmAddParam(prm,"nJeans", paramInt, &nJeans, sizeof(int),
"nJeans",
"<no. particle masses inside Jeans mass for good resolution> = 50");
dJeansConstant = 4/3.*pow(M_PI,2.5);
prmAddParam(prm,"dJeansConstant", paramDouble, &dJeansConstant,
sizeof(double), "dJeansConstant",
"<Constant to multiply c_s^3 G^-3/2 rho^-1/2 to get Jeans Mass > = 4/3 pi^5/2");
bSinkFormJeans = 0;
prmAddParam(prm,"bSinkFormJeans",paramBool,&bSinkFormJeans,
sizeof(int),
"sinkformjeans","enable/disable sinks = -sinkformjeans");
bSinkFormDivV = 0;
prmAddParam(prm,"bSinkFormDivV",paramBool,&bSinkFormDivV,
sizeof(int),
"sinkformdivv","enable/disable sinks = -sinkformdivv");
dSinkFormDivVCoeff = 0;
prmAddParam(prm,"dSinkFormDivVCoeff",paramDouble,
&dSinkFormDivVCoeff,sizeof(double),
"sinkformdivvcoeff","sinks form div v coeff = -sinkformdivvcoeff");
bSinkFormDivAcc = 0;
prmAddParam(prm,"bSinkFormDivAcc",paramBool,&bSinkFormDivAcc,
sizeof(int),
"sinkformdivacc","enable/disable sinks = -sinkformdivacc");
dSinkFormDivAccCoeff = 0;
prmAddParam(prm,"dSinkFormDivAccCoeff",paramDouble,
&dSinkFormDivAccCoeff,sizeof(double),
"sinkformdivvacccoeff","sinks form div acc coeff = -sinkformdivacccoeff");
bSinkFormDV = 0;
prmAddParam(prm,"bSinkFormDV",paramBool,&bSinkFormDV,sizeof(int),
"sinkformdv","enable/disable sinks = -sinkformdv");
bSinkFormPotMin = 0;
prmAddParam(prm,"bSinkFormPotMin",paramBool,&bSinkFormPotMin,
sizeof(int),
"sinkformpotmin","enable/disable sinks = -sinkformpotmin");
dSinkFormDensity = -1.0;
prmAddParam(prm,"dSinkFormDensity",paramDouble,
&dSinkFormDensity,sizeof(double),
"sinkformdensity","sinks density = -sinkformdensity");
dSinkTimeEligible = 0;
prmAddParam(prm,"dSinkTimeEligible",paramDouble,
&dSinkTimeEligible,sizeof(double),
"sinktimeligible","sink time eligible = -sinktimeeligible");
bSinkFormSimple = 0;
prmAddParam(prm,"bSinkFormSimple",paramBool,&bSinkFormSimple,
sizeof(int),
"sinkformsimple","enable/disable sinks = -sinkformjeans");
nSinkFormMin = 0;
prmAddParam(prm,"nSinkFormMin",paramInt,&nSinkFormMin,
sizeof(int), "sinkformmin","nmin to form a sink");
}
#include "physconst.h"
///
/// @brief check sink parameters
///
void Sinks::CheckParams(PRM prm, struct parameters ¶m)
{
if (!prmSpecified(prm, "nSinkFormMin")) nSinkFormMin = param.nSmooth;
if (bBHSink) {
assert (prmSpecified(prm, "dMsolUnit") &&
prmSpecified(prm, "dKpcUnit"));
/* For BH sinks -- default to negative tform as sink indicator */
if(!prmSpecified(prm, "dSinkMassMin"))
dSinkMassMin = FLT_MAX;
if(!param.bDoGas) {
ckerr << "Warning: BH sinks set without enabling SPH" << endl;
ckerr << "Enabling SPH" << endl;
param.bDoGas = 1;
}
bDoSinks = 1;
/* Units of inverse time -- code units */
dBHSinkEddFactor = GCGS*4*M_PI*MHYDR/
(SIGMAT*LIGHTSPEED*dBHSinkEddEff)*param.dSecUnit;
/* c^2 times efficiency factor (ergs per g) -- code units */
dBHSinkFeedbackFactor = dBHSinkFeedbackEff
*dBHSinkEddEff*LIGHTSPEED*LIGHTSPEED/param.dErgPerGmUnit;
dKmPerSecUnit = sqrt(GCGS*param.dMsolUnit*MSOLG
/(param.dKpcUnit*KPCCM))/1e5 ;
}
if (bDoSinks) {
double testDelta;
if (prmSpecified(prm,"iSinkRung")) {
/* Find associated timestep for iSinkRung */
int iRung;
if (iSinkRung > param.iMaxRung)
iSinkRung = param.iMaxRung;
testDelta = param.dDelta;
for ( iRung = 0; iRung < iSinkRung ; iRung++ )
testDelta *= 0.5;
}
else if (prmSpecified(prm,"dDeltaSink")) {
/* Find associate Rung for dDeltaSink */
/* NB: dDeltaSink is always in code units */
iSinkRung = 0;
for ( testDelta = param.dDelta; testDelta > dDeltaSink ;
testDelta *= 0.5 ) {
iSinkRung++;
if (iSinkRung >= param.iMaxRung)
CkAbort("Sink timestep too small");
}
}
else {
testDelta = param.dDelta;
iSinkRung = 0;
}
CkPrintf("dDeltaSink (set): %g, effectively: %g = %g yrs, iSinkRung: %i\n",
dDeltaSink, testDelta,
testDelta*param.dSecUnit/SECONDSPERYEAR, iSinkRung );
dDeltaSink = testDelta;
}
}
///
/// @brief Output black hole orbit information
///
/// Calls TreePiece::outputBlackHoles()
///
void Main::outputBlackHoles(double dTime)
{
FILE *fp;
double dOutTime;
double dvFac;
if (param.sinks.bBHSink) {
auto file_name = make_formatted_string("%s.orbit", param.achOutName);
char const* achOutFile = file_name.c_str();
fp = CmiFopen(achOutFile,"a");
CkAssert(fp != NULL);
fprintf(fp, "%ld %g\n", nSink, dTime);
long lFPos = ftell(fp);
CmiFclose(fp);
if(param.csm->bComove) {
dOutTime = csmTime2Exp(param.csm, dTime);
dvFac = 1.0/(dOutTime*dOutTime);
}
else {
dOutTime = dTime;
dvFac = 1.0;
}
treeProxy[0].outputBlackHoles(achOutFile, dvFac, lFPos,
CkCallbackResumeThread());
}
}
///
/// @brief processor specific routine to output black hole orbits.
///
/// Note this is not parallel.
///
void
TreePiece::outputBlackHoles(const std::string& pszFileName, double dvFac,
long lFPos, const CkCallback &cb)
{
FILE *fp;
int nout;
fp = CmiFopen(pszFileName.c_str(),"a");
CkAssert(fp != NULL);
int ret = fseek(fp, lFPos, SEEK_SET);
CkAssert(ret == 0);
/*
** Write Stuff!
*/
for(unsigned int i = 1; i <= myNumParticles; ++i) {
GravityParticle *p = &myParticles[i];
if(p->isStar() && p->fTimeForm() < 0.0)
fprintf(fp, "%li %g %g %g %g %g %g %g %g\n", p->iOrder,
p->mass, p->position.x, p->position.y, p->position.z,
dvFac*p->velocity.x, dvFac*p->velocity.y,
dvFac*p->velocity.z, p->potential);
}
lFPos = ftell(fp);
nout = CmiFclose(fp);
CkAssert(nout == 0);
if(thisIndex != (int) numTreePieces - 1)
pieces[thisIndex + 1].outputBlackHoles(pszFileName, dvFac, lFPos, cb);
else
cb.send();
}
///
/// @brief Initial identify sinks
///
void Main::SetSink()
{
if (param.sinks.bDoSinks) {
CkReductionMsg *msg;
treeProxy.SetSink(param.sinks.dSinkMassMin,
CkCallbackResumeThread((void*&)msg));
nSink = *(int *)msg->getData();
if (verbosity) CkPrintf("Identified %d sink particles\n", nSink);
delete msg;
}
}
///
/// @brief Initial identify sinks; processor specific method.
///
void TreePiece::SetSink(double dSinkMassMin, const CkCallback &cb)
{
int nSink = 0;
for(unsigned int i = 1; i <= myNumParticles; ++i) {
GravityParticle *p = &myParticles[i];
#ifdef STARSINK
if ((TYPETest(p,TYPE_STAR)))
#else
if (p->isStar() && (p->fTimeForm() < 0 || p->mass > dSinkMassMin))
#endif
{
TYPESet(p,TYPE_SINK);
nSink++;
}
}
contribute(sizeof(int), &nSink, CkReduction::sum_int, cb);
}
void TreePiece::SinkStep(int iCurrSinkRung, int iKickRung,
const CkCallback &cb)
{
for(unsigned int i = 1; i <= myNumParticles; ++i) {
GravityParticle *p = &myParticles[i];
if(p->rung >= iKickRung && TYPETest(p, TYPE_SINK|TYPE_SINKING)) {
if(iCurrSinkRung > p->rung) p->rung = iCurrSinkRung;
}
}
contribute(cb);
}
///
/// @brief Form sink particles; main routine
/// @param dTime Current time.
/// @param dDelta Current timestep.
/// @param iKickRung Rung being kicked.
///
/// Calls TreePiece::formSinks.
void Main::FormSinks(double dTime, double dDelta, int iKickRung)
/* For momentum conservation, only particles on iKickRung or higher may
contribute mass to a sink.
This call is timed for just after those particles have completed a full KDK
*/
{
int64_t nStarOld,nSinkAdd;
double sec,dsec;
if (!param.sinks.bSinkForm) return;
sec = CkWallTimer();
nStarOld = nTotalStar;
param.sinks.dSinkCurrentDelta = dDelta;
param.sinks.iSinkCurrentRung = iKickRung;
int bJeans = param.sinks.bSinkFormJeans;
int bSimple = param.sinks.bSinkFormSimple;
double dJConst2 = param.sinks.dJeansConstant/param.sinks.nJeans;
dJConst2 *= dJConst2;
double dDensityCut = param.sinks.dSinkFormDensity/param.dGmPerCcUnit;
int bDensity = (dDensityCut > 0 ? 1 : 0);
CkReductionMsg *msg;
treeProxy.formSinks(bJeans, dJConst2, bDensity, dDensityCut, dTime,
iKickRung, bSimple,
CkCallbackResumeThread((void*&)msg));
int nCandidates = *(int *)msg->getData();
delete msg;
nSinkAdd = 0;
if (nCandidates) {
if (param.sinks.bSinkFormSimple) {
nSinkAdd = nCandidates;
nSink += nSinkAdd;
}
else {
/* Note: Only gas particles are combined into sinks */
SinkFormTestSmoothParams pSFT(TYPE_GAS, iKickRung, param.sinks);
double dfBall2OverSoft2 = 4.0*param.dhMinOverSoft*param.dhMinOverSoft;
treeProxy.startSmooth(&pSFT, 1, param.nSmooth, dfBall2OverSoft2,
CkCallbackResumeThread());
SinkFormSmoothParams pSF(TYPE_GAS, iKickRung, param.csm, dTime,
param.sinks);
treeProxy.startSmooth(&pSF, 1, param.nSmooth, dfBall2OverSoft2,
CkCallbackResumeThread());
addDelParticles();
nSink += nSinkAdd = nTotalStar-nStarOld;
}
}
if (verbosity)
CkPrintf("Sinks Formation: %i candidates +%d sinks\n",nCandidates,
nSinkAdd);
dsec = CkWallTimer() - sec;
CkPrintf("Sink Form Done (+%d total) Calculated, Wallclock: %f secs\n\n",
nSinkAdd, dsec);
}
///
/// @brief Form sink particles: per processor routine
///
void TreePiece::formSinks(int bJeans, double dJConst2, int bDensity,
double dDensityCut, double dTime, int iKickRung,
int bSimple, const CkCallback &cb)
{
int i;
int n = myNumParticles;
double Jval, Jvalmin = FLT_MAX;
int nCandidates = 0;
for(i = 1; i <= n; ++i) {
GravityParticle *p = &myParticles[i];
if(p->isGas() && p->rung >= iKickRung) {
#ifdef SINKING
if (TYPETest( p, TYPE_SINKING )) continue;
#endif
/* Jeans Mass compared to nJeans particle masses */
Jval = dJConst2*(p->c()*p->c()*p->c()*p->c()*p->c()*p->c())
/(p->mass*p->mass*p->fDensity);
if (Jval < Jvalmin) Jvalmin = Jval;
if ((bJeans && Jval < 1) ||
(bDensity && p->fDensity >= dDensityCut)) {
nCandidates++;
if (bSimple) {
TYPESet(p, TYPE_SINK); /* Is now a sink! */
p->fMetals()= -dTime;
}
else {
TYPESet(p, TYPE_SMOOTHACTIVE); /* Is now a candidate */
}
}
}
}
contribute(cb);
}
int SinkFormTestSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_GAS));
}
#define fSinkRating(_a) ((_a)->curlv()[0] )
#define fSinkRating_Ex(_a) ((_a)->curlv[0] )
#define iOrderSink(_a) (*((int64_t *) (&(_a)->curlv()[1])))
#define iOrderSink_Ex(_a) (*((int64_t *) (&(_a)->curlv[1])))
void SinkFormTestSmoothParams::initSmoothCache(GravityParticle *p)
{
fSinkRating(p) = -FLT_MAX;
iOrderSink(p) = -1;
#ifdef SINKING
TYPEReset(p, TYPE_NEWSINKING);
#endif
}
void SinkFormTestSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
/* Particle p1 belongs to candidate stored in iOrderSink of p1 initially but
switch to p2's if that candidate is denser */
#ifdef SINKING
if (TYPETest(p1, TYPE_SINKING)) return;
#endif
if (fSinkRating_Ex(p2) > fSinkRating(p1)) {
fSinkRating(p1) = fSinkRating_Ex(p2);
iOrderSink(p1) = iOrderSink_Ex(p2);
}
}
void SinkFormTestSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
int i;
double dSinkRadius2 = s.dSinkRadius*s.dSinkRadius,r2;
GravityParticle *q;
/* Apply Bate tests in next phase
For now just decide which sink the particle belongs to:
prefer joining a denser (or deeper potential) sink candidate
Must be local extremum to be a candidate or abort
Need to prevent double counting particles into two sinks
*/
for (i=0;i<nSmooth;++i) {
r2 = nnList[i].fKey;
if (r2 > 0 && r2 <= dSinkRadius2) {
q = nnList[i].p;
if (TYPETest( q, TYPE_GAS ) && ((s.bSinkFormPotMin ? q->potential < p->potential : q->fDensity > p->fDensity))) {
/* Abort without grabbing any particles -- this isn't an extremum particle */
/* CkPrintf("Sink aborted %d %g: Denser Neighbour %d %g\n",p->iOrder,p->fDensity,q->iOrder,q->fDensity);*/
return;
}
}
}
fSinkRating(p) = (s.bSinkFormPotMin ? -p->potential : p->fDensity); /* rate yourself */
/* printf("Sink %d %g: looking...\n",p->iOrder,p->fDensity);*/
/* Identify all nbrs as to be eaten unless those nbrs already
belong to a "higher rated" sink
Thus sink quality is rated by density or potential depth
*/
for (i=0;i<nSmooth;++i) {
r2 = nnList[i].fKey;
if (r2 > 0 && r2 <= dSinkRadius2) {
q = nnList[i].p;
#ifdef SINKING
if (TYPETest(q, TYPE_SINKING)) continue;
#endif
if (TYPETest( q, TYPE_GAS ) && q->rung >= s.iSinkCurrentRung) {
if (s.bSinkFormPotMin) {
if (-p->potential > fSinkRating(q)) {
fSinkRating(q) = -p->potential;
iOrderSink(q) = p->iOrder; /* Particle q belongs to sink p */
}
}
else {
if (p->fDensity > fSinkRating(q)) {
fSinkRating(q) = p->fDensity;
iOrderSink(q) = p->iOrder; /* Particle q belongs to sink p */
}
}
}
}
}
}
void tred3(double a[][3], double d[], double e[])
{
CkAbort("Matrix diagnolization not implemented");
}
void tqli3(double d[], double e[], double z[][3])
{
CkAbort("Matrix diagnolization not implemented");
}
int SinkFormSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_GAS));
}
void SinkFormSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if (!(TYPETest( p1, TYPE_DELETED )) && TYPETest( p2, TYPE_DELETED ) ) {
p1->mass = p2->mass;
deleteParticle(p1);
}
}
void SinkFormSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
int i,j,nEaten,nEatNow;
double mtot,im,Ek,Eth,Eg,dvx,dv2,vsink[3];
double diva,ih2,r2,rs1;
GravityParticle *q,*q1,*q2;
GravityParticle sinkp;
/* You are accreted */
if ( iOrderSink(p) != -1 ) {
/* What if it was temporarily a candidate -- are it's nbrs confused?
Shouldn't be -- they should go to a new better candidate if they are near one */
return;
}
#ifdef SINKING
assert(!TYPETest(p,TYPE_SINKING));
#endif
iOrderSink(p) = p->iOrder;
nEaten=0;
mtot = 0;
for (i=0;i<nSmooth;++i) {
q1 = nnList[i].p;
if (iOrderSink(q1) == p->iOrder) {
nEaten++;
mtot += q1->mass;
}
}
if (nEaten < s.nSinkFormMin) {
return;
}
if (s.bSinkFormDivAcc) {
diva = 0;
ih2 = 1./(s.dSinkRadius*s.dSinkRadius);
for (i=0;i<nSmooth;++i) {
q1 = nnList[i].p;
if (iOrderSink(q1) == p->iOrder) {
double dax,day,daz,dx,dy,dz;
r2 = nnList[i].fKey*ih2; /* Calculate div a (un-normalized) */
if (r2 < 4) {
rs1 = DKERNEL(r2);
rs1 *= q1->mass;
dx = nnList[i].dx.x;
dy = nnList[i].dx.y;
dz = nnList[i].dx.z;
dax = (-p->treeAcceleration[0]
+ q1->treeAcceleration[0]);
day = (-p->treeAcceleration[1]
+ q1->treeAcceleration[1]);
daz = (-p->treeAcceleration[2]
+ q1->treeAcceleration[2]);
diva += (dax*dx+day*dy+daz*dz)*rs1;
}
}
}
if (diva > 0) return;
}
if (s.bSinkFormDV || s.bSinkFormDivV) {
double dvdx[3][3];
double d[3],e[3];
double vFac = 1./(a*a); /* converts v to xdot */
dvdx[0][0]=0; dvdx[0][1]=0; dvdx[0][2]=0;
dvdx[1][0]=0; dvdx[1][1]=0; dvdx[1][2]=0;
dvdx[2][0]=0; dvdx[2][1]=0; dvdx[2][2]=0;
ih2 = 1./(s.dSinkRadius*s.dSinkRadius);
for (i=0;i<nSmooth;++i) {
q1 = nnList[i].p;
if (iOrderSink(q1) == p->iOrder) {
double dvx,dvy,dvz,dx,dy,dz;
r2 = nnList[i].fKey*ih2; /* Calculate dv_i/dx_j (un-normalized) */
if (r2 < 4) {
rs1 = DKERNEL(r2);
rs1 *= q1->mass;
dx = nnList[i].dx.x;
dy = nnList[i].dx.y;
dz = nnList[i].dx.z;
dvx = ((-p->vPred()[0] + q1->vPred()[0])*vFac - dx*H)*rs1; /* NB: dx = px - qx */
dvy = ((-p->vPred()[1] + q1->vPred()[1])*vFac - dy*H)*rs1;
dvz = ((-p->vPred()[2] + q1->vPred()[2])*vFac - dz*H)*rs1;
dvdx[0][0] += dvx*dx;
dvdx[0][1] += dvx*dy;
dvdx[0][2] += dvx*dz;
dvdx[1][0] += dvy*dx;
dvdx[1][1] += dvy*dy;
dvdx[1][2] += dvy*dz;
dvdx[2][0] += dvz*dx;
dvdx[2][1] += dvz*dy;
dvdx[2][2] += dvz*dz;
}
}
}
dvdx[0][1] = 0.5*(dvdx[0][1]+dvdx[1][0]); dvdx[1][0]=dvdx[0][1];
dvdx[2][1] = 0.5*(dvdx[2][1]+dvdx[1][2]); dvdx[1][2]=dvdx[2][1];
dvdx[2][0] = 0.5*(dvdx[2][0]+dvdx[0][2]); dvdx[0][2]=dvdx[2][0];
if (s.bSinkFormDivV) {
if (dvdx[0][0]+dvdx[1][1]+dvdx[2][2] > 0) return; /* not converging */
}
if (s.bSinkFormDV) {
tred3(dvdx,d,e);
tqli3(d,e,dvdx);
if (d[0] > 0 || d[1] > 0 || d[2] > 0) return; /* not converging in all directions */
}
}
Ek = 0;
Eth = 0;
Eg = 0;
vsink[0] = 0; vsink[1] = 0; vsink[2] = 0;
for (i=0;i<nSmooth;++i) {
q1 = nnList[i].p;
if (iOrderSink(q1) == p->iOrder) {
dvx = q1->velocity[0]-p->velocity[0];
vsink[0] += q1->mass*dvx;
dv2 = dvx*dvx;
dvx = q1->velocity[1]-p->velocity[1];
vsink[1] += q1->mass*dvx;
dv2 += dvx*dvx;
dvx = q1->velocity[2]-p->velocity[2];
vsink[2] += q1->mass*dvx;
dv2 += dvx*dvx;
Ek += 0.5*q1->mass*dv2;
Eth += q1->mass*q1->u();
for (j=i+1;j<nSmooth;j++) {
q2 = nnList[j].p;
if (iOrderSink(q2) == p->iOrder) {
dvx = q1->position[0]-q2->position[0];
dv2 = dvx*dvx;
dvx = q1->position[1]-q2->position[1];
dv2 += dvx*dvx;
dvx = q1->position[2]-q2->position[2];
dv2 += dvx*dvx;
Eg -= q1->mass*q2->mass/sqrt(dv2);
}
}
}
}
Ek -= 0.5*(vsink[0]*vsink[0]+vsink[1]*vsink[1]+vsink[2]*vsink[2])/mtot;
/* Apply Bate tests here --
1. thermal energy < 1/2 |Grav|,
2. thermal + rot E < |Grav|,
3. total E < 0 (implies 2.)
4. div.acc < 0 (related to rate of change of total E I guess) (optional see above)
5. NEW eigenvalues of dvi/dxj all negative (Federrath) (optional see above)
6. NEW div.v < 0 (Whitworth group)
*/
if (Eth < 0.5*fabs(Eg) && Ek + Eth + Eg < 0) {
/* Sink approved */
if (s.dSinkTimeEligible > 0) {
/* Check that formation conditions persisted for minimum time */
/* This cannot work in its simplest form because the sink candidate
for a given density peak changes randomly from step to step as
densities etc.. fluctuate. Need to code in some way to continuously
refer to the same pre-sink blob (inherit properties of particles inside 1/2 SinkRadius?) */
if (p->fTimeForm() > 1e36) {
p->fTimeForm() = dTime;
return;
}
if (dTime-p->fTimeForm() < s.dSinkTimeEligible) return;
}
GravityParticle sinkp = *p;
sinkp.position[0] = 0;
sinkp.position[1] = 0;
sinkp.position[2] = 0;
sinkp.velocity[0] = 0;
sinkp.velocity[1] = 0;
sinkp.velocity[2] = 0;
sinkp.treeAcceleration[0] = 0;
sinkp.treeAcceleration[1] = 0;
sinkp.treeAcceleration[2] = 0;
sinkp.u() = 0;
#ifdef STARSINK
SINK_Lx(&sinkp) = 0;
SINK_Ly(&sinkp) = 0;
SINK_Lz(&sinkp) = 0;
#endif
sinkp.mass = 0;
#ifdef SINKING
sinkp.fTrueMass = 0;
sinkp.iSinkingOnto = 0; /* Must start at zero -- it records the number of assoc. sinking particles */
#endif
nEatNow = 0;
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
r2 = nnList[i].fKey;
if (iOrderSink(q) == p->iOrder) {
#ifdef SINKING
if (r2 > smf->dSinkMustAccreteRadius*smf->dSinkMustAccreteRadius) {
/* Make sinking -- can't force sink choice b/c sink iOrder isn't known yet */
TYPESet(q, TYPE_NEWSINKING);
}
else
#endif /*SINKING*/
{
double mq,rx,ry,rz,vx,vy,vz;
mq = q->mass;
rx = q->position[0]; ry = q->position[1];
rz = q->position[2];
vx = q->velocity[0]; vy = q->velocity[1];
vz = q->velocity[2];
nEatNow++;
sinkp.position[0] += mq*rx;
sinkp.position[1] += mq*ry;
sinkp.position[2] += mq*rz;
sinkp.velocity[0] += mq*vx;
sinkp.velocity[1] += mq*vy;
sinkp.velocity[2] += mq*vz;
sinkp.treeAcceleration[0] += mq*q->treeAcceleration[0];
sinkp.treeAcceleration[1] += mq*q->treeAcceleration[1];
sinkp.treeAcceleration[2] += mq*q->treeAcceleration[2];
sinkp.u() += q->mass*q->u();
#ifdef STARSINK
SINK_Lx(&sinkp) += mq*(ry*vz - rz*vy); /* Add to total L */
SINK_Ly(&sinkp) += mq*(rz*vx - rx*vz);
SINK_Lz(&sinkp) += mq*(rx*vy - ry*vx);
#endif
sinkp.mass += mq;
#ifdef SINKING
sinkp.fTrueMass += mq;
#endif
if (p!=q) {
q->mass = 0;
deleteParticle(q);
}
}
}
}
im = 1/sinkp.mass;
sinkp.position[0] *= im;
sinkp.position[1] *= im;
sinkp.position[2] *= im;
sinkp.velocity[0] *= im;
sinkp.velocity[1] *= im;
sinkp.velocity[2] *= im;
sinkp.treeAcceleration[0] *= im;
sinkp.treeAcceleration[1] *= im;
sinkp.treeAcceleration[2] *= im;
sinkp.u() *= im;
#ifdef STARSINK
/* Store Internal angular momentum only */
SINK_Lx(&sinkp) -= sinkp.fMass*(sinkp.r[1]*sinkp.v[2] - sinkp.r[2]*sinkp.v[1]);
SINK_Ly(&sinkp) -= sinkp.fMass*(sinkp.r[2]*sinkp.v[0] - sinkp.r[0]*sinkp.v[2]);
SINK_Lz(&sinkp) -= sinkp.fMass*(sinkp.r[0]*sinkp.v[1] - sinkp.r[1]*sinkp.v[0]);
#endif
TYPEClear(&sinkp);
TYPESet(&sinkp,TYPE_SINK|TYPE_STAR);
sinkp.fTimeForm() = -dTime; /* -ve time is sink indicator */
#ifdef SINKEXTRADATA
for(j = 0; j < 3; j++) {
sinkp.rForm[j] = sinkp.r[j];
sinkp.vForm[j] = sinkp.v[j];
}
#endif
printf("Sink Formed %ld %g: np (%d) %d Mass (%g) %g Ek %g Eth %g Eg %g, %g %g\n",p->iOrder,p->fDensity,nEaten,nEatNow,mtot,sinkp.mass,Ek,Eth,Eg,4/3.*pow(M_PI,2.5)/50.*sqrt((p->c()*p->c()*p->c()*p->c()*p->c()*p->c())/(p->mass*p->mass*p->fDensity)),pow(Eth/fabs(Eg),1.5) );
#ifndef SINKING
assert(fabs(sinkp.mass/mtot-1) < 1e-4);
#else
/* printf("Sink Made %g %g\n",sinkp.fMass,sinkp.fTrueMass);*/
assert(fabs(sinkp.fMass/sinkp.fTrueMass-1) < 1e-7);
assert(!TYPETest(p, TYPE_SINKING));
#endif
p->mass = 0;
deleteParticle(p);
tp->newParticle(&sinkp);
}
else {
CkPrintf("Sink Failed Tests %ld %g: np (%d) Mass (%g) Ek %g Eth %g Eg %g, %g %g\n",p->iOrder,p->fDensity,nEaten,mtot,Ek,Eth,Eg, 4/3.*pow(M_PI,2.5)/50.*sqrt((p->c()*p->c()*p->c()*p->c()*p->c()*p->c())/(p->mass*p->mass*p->fDensity)),pow(Eth/fabs(Eg),1.5) );
}
}
/// Process sink particles
void Main::doSinks(double dTime, double dDelta, int iKickRung)
/* For momentum conservation, only particles on iKickRung or higher may
contribute mass to a sink
This call is timed for just after those particles have completed a full KDK
*/
{
double sec,sec1,dsec;
int nAccreted;
/* I assume sink creation is rarer so the tree will be ok after this call most of the time */
FormSinks(dTime, dDelta, iKickRung );
if (param.sinks.bDoSinks == 0) return;
if (nSink == 0) return;
if (param.sinks.bBHSink && dDelta <= 0.0) return;
sec = CkWallTimer();
nAccreted = nTotalSPH;
if(nActiveGrav > 0) {
param.sinks.dSinkCurrentDelta = dDelta;
param.sinks.iSinkCurrentRung = iKickRung;
double dfBall2OverSoft2 = 4.0*param.dhMinOverSoft*param.dhMinOverSoft;
if (param.sinks.bBHSink) {
/* Smooth Bondi-Hoyle Accretion: radius set by nSmooth */
BHDensitySmoothParams pBHDen(TYPE_GAS, iKickRung, param.csm,
dTime, param.dDelta, param.sinks);
treeProxy.startSmooth(&pBHDen, 1, param.nSmooth, dfBall2OverSoft2,
CkCallbackResumeThread());
BHAccreteSmoothParams pBHAcc(TYPE_GAS, iKickRung, param.csm,
dTime, param.dDelta, param.sinks,
param.dConstGamma,
param.stfm->dMinGasMass);
treeProxy.startReSmooth(&pBHAcc, CkCallbackResumeThread());
/* Search for BHs that need merging */
CkReductionMsg *msgCnt;
treeProxy.countType(TYPE_SINK,
CkCallbackResumeThread((void *&)msgCnt));
nSink = *(int *) msgCnt->getData();
delete msgCnt;
if(nSink > 0) {
if(verbosity)
CkPrintf("Distribute Deleted gas\n");
DistDeletedGasSmoothParams pDGas(TYPE_GAS, 0);
treeProxy.startReSmooth(&pDGas, CkCallbackResumeThread());
}
if (nSink > 1) { /* no need to merge if there is only one! */
BHIdentifySmoothParams pBHID(TYPE_SINK, iKickRung, param.csm,
dTime, param.sinks);
int nSmooth = nSink < 4 ? nSink : 4;
treeProxy.startSmooth(&pBHID, 1, nSmooth, dfBall2OverSoft2,
CkCallbackResumeThread());
BHSinkMergeSmoothParams pBHM(TYPE_SINK, iKickRung, param.csm,
dTime, param.sinks,
param.dConstGamma);
treeProxy.startReSmooth(&pBHM, CkCallbackResumeThread());
}
}
else {
/* Fixed Radius Accretion: particle by particle (cf. Bate) */
SinkAccreteTestSmoothParams pSinkTest(TYPE_GAS, iKickRung, dTime,
param.sinks);
treeProxy.startSmooth(&pSinkTest, 1, param.nSmooth,
dfBall2OverSoft2, CkCallbackResumeThread());
#ifdef SINKINGAVERAGE
msrActiveType(msr,TYPE_NEWSINKING, TYPE_SMOOTHACTIVE);
msrSmooth(msr, dTime, SMX_SINKINGAVERAGE,0);
SinkingAverageSmoothParams pSinkAvg(TYPE_GAS, iKickRung, sinks);
treeProxy.startSmooth(&pSinkAvg, 1, param.nSmooth,
dfBall2OverSoft2, CkCallbackResumeThread());
#endif
SinkAccreteSmoothParams pSinkAcc(TYPE_GAS, iKickRung, dTime,
param.sinks);
treeProxy.startSmooth(&pSinkAcc, 1, param.nSmooth,
dfBall2OverSoft2, CkCallbackResumeThread());
}
addDelParticles();
if (param.sinks.bBHSink) { /* reset nSinks in case of merger */
CkReductionMsg *msgCnt;
treeProxy.countType(TYPE_SINK,
CkCallbackResumeThread((void *&)msgCnt));
nSink = *(int *) msgCnt->getData();
delete msgCnt;
}
}
nAccreted -= nTotalSPH;
sec1 = CkWallTimer();
dsec = sec1 - sec;
CkPrintf("Sinks Done (%d accreted) Calculated, Wallclock: %f secs\n\n",
nAccreted,dsec);
}
int BHDensitySmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
void BHDensitySmoothParams::initTreeParticle(GravityParticle *p)
{
p->curlv()[1] = 0.0; /* total mass change */
p->curlv()[0] = p->mass; /* initial mass */
p->curlv()[2] = 0.0; /* flag if accreted */
}
/* Cached Tree Active particles */
void BHDensitySmoothParams::initSmoothCache(GravityParticle *p)
{
/*
* Original particle curlv's is trashed (JW)
*/
p->curlv()[1] = 0.0; /* total mass change */
p->curlv()[0] = p->mass; /* initial mass */
p->curlv()[2] = 0.0; /* flag if accreted */
}
void BHDensitySmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if (p2->curlv[2] > p1->curlv()[2]) { /* only important if particle has been selected JMB 6/24/09
* but necessary - if particle has 2 BH neighbors it can be
* overwritten */
p1->curlv()[0] = p2->curlv[0];
p1->curlv()[1] = p2->curlv[1]; /* total mass change */
p1->curlv()[2] = p2->curlv[2];
}
/* this way the most recent one wins JMB 5/12/09 */
/* don't sum up curlv[1], only let one BH eating be recorded.
curlv[1] will never be greater than curlv[0] this way */
}
/*
* Calculate parameters for BH accretion for use in the BHSinkAccrete
* function. This is done separately to even competition between
* neighboring Black Holes.
* Gas particles are ID'd here to be eaten in BHSinkAccrete.
*/
void BHDensitySmoothParams::fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nnList)
{
GravityParticle *q = NULL;
double ih2,r2,rs,fDensity;
double v[3],cs,fW,dv2,dv;
double mdot, mdotEdd, mdotCurr, dmq, dtEff;
int i, naccreted, ieat, ivmin;
double mdotsum, weat;
double weight,wrs;
double aFac, dCosmoDenFac,dCosmoVel2Fac,aDot;
double dvmin,dvx,dvy,dvz,dvcosmo; /* measure mindeltav */
assert(p->rung >= s.iSinkCurrentRung);
p->dDeltaM() = 0.0;
naccreted = 0;
aFac = a;
dCosmoDenFac = aFac*aFac*aFac;
dCosmoVel2Fac = aFac*aFac;
aDot = aFac*H; /* da/dt */
/* CosmoVel2Fac converts cannonical velocity^2 to physical velocity^2
Cosmo sims tend to use cannonical velocities
v_phys = v_cannonical/a + adot*r
JMB 9/1/12*/
mdotsum = 0.;
weat = -1e37;
ih2 = invH2(p);
fDensity = 0.0; cs = 0;
v[0] = 0; v[1] = 0; v[2] = 0;
/* for mindv accretion method JMB 8/5/10 */
if(s.bBHMindv == 1){
dvmin = FLT_MAX;
for (i=0;i<nSmooth;++i) {
dvx = (-p->velocity[0] + nnList[i].p->vPred()[0])/aFac - aDot*nnList[i].dx.x;
dvy = (-p->velocity[1] + nnList[i].p->vPred()[1])/aFac - aDot*nnList[i].dx.y;
dvz = (-p->velocity[2] + nnList[i].p->vPred()[2])/aFac - aDot*nnList[i].dx.z;
dvcosmo = sqrt(dvx*dvx + dvy*dvy + dvz*dvz);
if (dvcosmo < dvmin) {
dvmin=dvcosmo;
ivmin = i;
}
}
}
for (i=0;i<nSmooth;++i) {
double fDist2 = nnList[i].fKey;
r2 = fDist2*ih2;
rs = KERNEL(r2, nSmooth);
q = nnList[i].p;
q->curlv()[0] = q->mass;
q->curlv()[1] = 0;
q->curlv()[2] = 0;
/* Sink should NOT be part of this list */
assert(TYPETest(q,TYPE_GAS));
fW = rs*q->mass;
if(s.bBHMindv == 1) weight = rs*pow(q->c()*q->c()+(dvmin*dvmin),-1.5)/dCosmoDenFac;
else {
dvx = (-p->velocity[0]+q->vPred()[0])/aFac - aDot*nnList[i].dx.x;
dvy = (-p->velocity[1]+q->vPred()[1])/aFac - aDot*nnList[i].dx.y;
dvz = (-p->velocity[2]+q->vPred()[2])/aFac - aDot*nnList[i].dx.z;
dvcosmo = sqrt(dvx*dvx+dvy*dvy+dvz*dvz);
weight = rs*pow(q->c()*q->c()+dvcosmo*dvcosmo,-1.5)/dCosmoDenFac; /* weight particles by mdot quantities */
/* cosmo factors put in 7/7/09 JMB */
}
if (weight > weat) {
weat = weight;
ieat = i;
wrs = rs;
}
mdotsum += weight*q->mass;
fDensity += fW;
v[0] += fW*q->velocity[0];
v[1] += fW*q->velocity[1];
v[2] += fW*q->velocity[2];
cs += fW*q->c();
}
dv2 = 0;
for (i=0;i<3;i++) {
dv = v[i]/fDensity-p->velocity[i];
dv2 += dv*dv;
}
dv2 /= dCosmoVel2Fac;
/*
* Store results in particle.
*/
cs = cs/fDensity;
p->fDensity = fDensity = M_1_PI*sqrt(ih2)*ih2*fDensity;
fDensity /= dCosmoDenFac;
if(s.bBHMindv == 1)
CkPrintf("BHSink %ld: Time: %g Selected Particle Density: %g SPH: %g C_s: %g dv: %g dist: %g\n",
p->iOrder,dTime,nnList[ieat].p->fDensity, wrs,
nnList[ieat].p->c(),dvmin, sqrt(nnList[ieat].fKey));
else CkPrintf("BHSink %ld: Time: %g Selected Particle Density: %g SPH: %g C_s: %g dv: %g dist: %g\n",
p->iOrder,dTime,nnList[ieat].p->fDensity,wrs,
nnList[ieat].p->c(),sqrt(dv2),sqrt(nnList[ieat].fKey));
CkPrintf("BHSink %ld: Time: %g Mean Particle Density: %g C_s: %g dv: %g\n",p->iOrder,dTime,fDensity,cs,sqrt(dv2));
mdot = mdotsum*s.dBHSinkAlpha*4*M_PI*p->mass*p->mass*M_1_PI*sqrt(ih2)*ih2; /* new mdot! */
/* Eddington Limit Rate */
mdotEdd = s.dBHSinkEddFactor*p->mass;
CkPrintf("BHSink %ld: Time: %.8f mdot (BH): %g mdot (Edd): %g a: %g\n",
p->iOrder,dTime,mdot,mdotEdd, a);
if (mdot > mdotEdd) mdot = mdotEdd;
mdotCurr = p->dMDot() = mdot;
weight = -1e37;
weat = -1e37;
for (;;) {
q = NULL;
for (i=0;i<nSmooth;++i) {
double fDist2 = nnList[i].fKey;
r2 = fDist2;
if(TYPETest(nnList[i].p,TYPE_DELETED)) continue;
/* don't accrete a deleted particle! JMB 10/1/08 */
rs = KERNEL(r2, nSmooth);
fW = rs*nnList[i].p->mass;
double fBall = nnList[i].p->fBall;
if(!s.bBHAccreteAll && r2 > 0.25*fBall*fBall) continue;
/* don't accrete gas that doesn't have the BH
* in its smoothing length JMB 10/22/08 */
/* make this an optional parameter JMB 9/21/12 */
/* if (nnList[i].p->rung < s.iSinkCurrentRung) continue; /* JMB 7/9/09 */
if(s.bBHMindv == 1) weight = rs*pow(nnList[i].p->c()*nnList[i].p->c()+(dvmin*dvmin),-1.5)/dCosmoDenFac;
else {
dvx = (-p->velocity[0]+nnList[i].p->vPred()[0])/aFac - nnList[i].dx.x*aDot;
dvy = (-p->velocity[1]+nnList[i].p->vPred()[1])/aFac - nnList[i].dx.y*aDot;
dvz = (-p->velocity[2]+nnList[i].p->vPred()[2])/aFac - nnList[i].dx.z*aDot;
dvcosmo = sqrt(dvx*dvx+dvy*dvy+dvz*dvz);
weight = rs*pow(nnList[i].p->c()*nnList[i].p->c()+dvcosmo*dvcosmo,-1.5)/dCosmoDenFac; /* weight particles by mdot quantities */
/* cosmo factors put in 7/7/09 JMB */
}
if (weight > weat && nnList[i].p->curlv()[2] != p->iOrder) {
weat = weight;
q = nnList[i].p;
}
}
assert(q != p); /* shouldn't happen because p isn't in
tree */
if(q != NULL) {
#if 0
/* Timestep for accretion is larger of sink and victim timestep */
iRung = q->rung;
if (iRung > p->rung) iRung = p->rung;
#endif
// dtEff = s.dSinkCurrentDelta;
dtEff = RungToDt(dDelta, p->rung);
/* If victim has unclosed kick -- don't actually take the mass
If sink has unclosed kick we shouldn't even be here!
When victim is active use his timestep if longer
Statistically expect to get right effective mdot on average */
/* JMB 9/18/12 -- Actually, the method described above
does NOT get the right mdot on average, it
underestimates mdot. Instead I have removed all
timestep criteria from BH accretion. Thus, momentum
is no longer strictly conserved. HOWEVER, the
amounts of momentum involved are tiny, likely
comparable to the errors intrinsic to the code.
Thus, we shall not stress about it. */
dmq = mdotCurr*dtEff;
}
else {
if(naccreted != 0) CkPrintf("BHSink %ld: only %i accreted, %g uneaten mass\n",p->iOrder,naccreted,mdotCurr*dtEff);
/* JMB 6/29/09 */
if(naccreted == 0) p->dDeltaM() = 0.0; /* should be anyway */
return;
}
assert(q != NULL);
/* We have our victim */
if (dmq < q->mass) {
mdotCurr = 0.0;
}
else {
mdotCurr -= mdotCurr*(q->mass/dmq); /* need an additional victim */
assert(mdotCurr >= 0.0);
dmq = q->mass;
}
q->curlv()[2] = p->iOrder; /* flag for pre-used victim particles */
/* temporarily store mass lost -- to later check for double dipping */
q->curlv()[1] += dmq;
p->dDeltaM() += dmq;
naccreted += 1;
weat = -1e37; /* reset so other particles can be selected JMB 7/9/09 */
CkPrintf("BHSink %ld: Time %g %ld dmq %g %g %g\n",p->iOrder,dTime,q->iOrder,dmq,q->curlv()[1],p->dDeltaM());
if (mdotCurr == 0.0) break;
}
}
int BHAccreteSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
void BHAccreteSmoothParams::initTreeParticle(GravityParticle *p1)
{
/* Convert energy and metals to non-specific quantities (not per mass)
* to make it easier to divvy up BH energy.
*/
if(TYPETest(p1, TYPE_GAS)) p1->u() *= p1->mass;
if(TYPETest(p1, TYPE_GAS)) p1->uPred() *= p1->mass;
}
/* Cached Tree Active particles */
void BHAccreteSmoothParams::initSmoothCache(GravityParticle *p)
{
if (TYPETest(p, iType)) p->u() = 0.0; /*added 6/10/08*/
if (TYPETest(p, iType)) p->uPred() = 0.0;
}
static inline double max(double x, double y)
{
if (x > y) return x;
else return y;
}
void BHAccreteSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if (!(TYPETest( p1, TYPE_DELETED )) &&
TYPETest( p2, TYPE_DELETED ) ) {
p1->mass = p2->mass;
deleteParticle(p1);
}
else if (TYPETest(p1, iType)) { // TreeActive
/*
* See kludgery notice above: record eaten mass in original
* particle.
*/
double fEatenMass = p2->curlv[0] - p2->mass;
p1->mass -= fEatenMass;
if(p1->mass < p1->curlv()[0]*1e-3) {
/* This could happen if BHs on two
different processors are eating
gas from a third processor */
CkError("ERROR: Overeaten gas particle %ld: %g %g\n",
p1->iOrder, p1->mass, fEatenMass);
if (!(TYPETest( p1, TYPE_DELETED ))) {
deleteParticle( p1 );
}
return;
}
if(!(TYPETest(p1,TYPE_DELETED))) { /*added 8/21/08 to check -0 mass particles */
CkAssert(p1->mass > 0.0);
}
p1->u() += p2->u; /*added 6/10/08 for BH blastwave FB*/
p1->uPred() += p2->uPred;
p1->fTimeCoolIsOffUntil() = max( p1->fTimeCoolIsOffUntil(),
p2->fTimeCoolIsOffUntil );
/* propagate FB time JMB 2/24/10 */
p1->dTimeFB() = max(p1->dTimeFB(), p2->dTimeFB);
}
}
void BHAccreteSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
GravityParticle *q = NULL;
double ih2,r2,rs;
double mdot, mdotCurr, dmAvg, dm, dmq, dE, ifMass, dtEff;
double fNorm,fNorm_new,f2h2;
int i, counter,naccreted,ieat;
double weat;
double weight,fbweight; /* weight is for accretion, fbweight is for FB */
double aFac, dCosmoDenFac,dCosmoVel2Fac,aDot;
double dvmin, dvx,dvy,dvz,dvcosmo;
int ivmin;
weat = -1e37;
aFac = a;
dCosmoDenFac = aFac*aFac*aFac;
dCosmoVel2Fac = aFac*aFac;
aDot = aFac*H;
mdot = p->dMDot();
#if 0
if (p->dDeltaM() == 0.0) {
dtEff = s.dSinkCurrentDelta*pow(0.5,p->rung-s.iSinkCurrentRung);
dmAvg = mdot*dtEff;
CkPrintf("BHSink %ld: Delta: %g dm: 0 ( %g ) (victims on wrong step)\n",p->iOrder,dtEff,dmAvg);
return;
}
#endif
mdotCurr = mdot;
dm = 0;
naccreted = 0;
ih2 = invH2(p);
for (i=0;i<nSmooth;++i){
double r2min = FLT_MAX;
q = NULL;
if (nnList[i].p->curlv()[2] == 0) continue;
if(TYPETest(nnList[i].p,TYPE_DELETED)) continue;
if (nnList[i].p->curlv()[2] == p->iOrder) q = nnList[i].p;
else continue; /* won't work if selected by another BH */
assert(q != NULL);
r2min = nnList[i].fKey;
#if 0
/* Timestep for accretion is larger of sink and victim timestep */
iRung = q->rung;
if (iRung > p->rung) iRung = p->rung;
#endif
// dtEff = s.dSinkCurrentDelta;
dtEff = RungToDt(dDelta, p->rung);
/* If victim has unclosed kick -- don't actually take the mass
If sink has unclosed kick we shouldn't even be here!
When victim is active use his timestep if longer
Statistically expect to get right effective mdot on average */
/* JMB 9/18/12 - see comment in BHSinkDensity above for
why the previous comment is not actually what is done. */
dmq = mdotCurr*dtEff;
/* update mdotCurr */
if( fabs(dmq - q->curlv()[1])/dmq < 0.0001) mdotCurr = 0.0;
/* some floating point issues here JMB 7/9/09 */
else {
mdotCurr -= q->curlv()[1]/dtEff;
dmq = q->curlv()[1];
}
CkAssert(mdotCurr >= 0.0);
CkPrintf("BHSink %ld: Time %g %ld dmq %g %g %g\n",p->iOrder, dTime,
q->iOrder,dmq,q->curlv()[1],p->dDeltaM());
ifMass = 1./(p->mass + dmq);
/* to record angular momentum JMB 11/9/10 */
CkPrintf("BHSink %ld: Gas: %ld dx: %g dy: %g dz: %g \n",
p->iOrder,q->iOrder,p->position[0]-q->position[0],
p->position[1]-q->position[1],
p->position[2]-q->position[2]);
CkPrintf("BHSink %ld: Gas: %ld dvx: %g dvy: %g dvz: %g \n",
p->iOrder,q->iOrder,p->velocity[0]-q->velocity[0],
p->velocity[1]-q->velocity[1],
p->velocity[2]-q->velocity[2]);
/* Adjust sink properties (conserving momentum etc...) */
p->position[0] = ifMass*(p->mass*p->position[0]
+dmq*q->position[0]);
p->position[1] = ifMass*(p->mass*p->position[1]
+dmq*q->position[1]);
p->position[2] = ifMass*(p->mass*p->position[2]
+dmq*q->position[2]);
p->velocity[0] = ifMass*(p->mass*p->velocity[0]
+dmq*q->velocity[0]);
p->velocity[1] = ifMass*(p->mass*p->velocity[1]
+dmq*q->velocity[1]);
p->velocity[2] = ifMass*(p->mass*p->velocity[2]
+dmq*q->velocity[2]);
p->treeAcceleration[0] = ifMass*(p->mass*p->treeAcceleration[0]
+dmq*q->treeAcceleration[0]);
p->treeAcceleration[1] = ifMass*(p->mass*p->treeAcceleration[1]
+dmq*q->treeAcceleration[1]);
p->treeAcceleration[2] = ifMass*(p->mass*p->treeAcceleration[2]
+dmq*q->treeAcceleration[2]);
p->fStarMetals() = ifMass*(p->mass*p->fStarMetals()+dmq*q->fMetals());
p->mass += dmq;
dm += dmq;
q->mass -= dmq;
/*assert(q->mass >= 0.0);*/
naccreted += 1; /* running tally of how many are accreted JMB 10/23/08 */
CkPrintf("BHSink %ld: Time %g dist2: %ld %g gas smooth: %g eatenmass %g \n",
p->iOrder,dTime,q->iOrder,r2min,q->fBall*q->fBall,dmq);
if (q->mass <= dMinGasMass) {
if(!(TYPETest(q,TYPE_DELETED))) deleteParticle(q);
/* Particles are getting deleted twice, which is messing
up the bookkeeping. I think this will make them only be
deleted once. JMB 9/23/08*/
}
}
if (mdotCurr == 0.0) goto dofeedback;
else{
/* go looking for more particles, sharing must happen */
/****** SHARING CODE HERE ***********/
if(s.bBHMindv == 1){
dvmin = FLT_MAX;
for (i=0;i<nSmooth;++i) {
dvx = (-p->velocity[0] + nnList[i].p->vPred()[0])-aDot*nnList[i].dx.x;
dvy = (-p->velocity[1] + nnList[i].p->vPred()[1])-aDot*nnList[i].dx.y;
dvz = (-p->velocity[2] + nnList[i].p->vPred()[2])-aDot*nnList[i].dx.z;
dvcosmo = sqrt(dvx*dvx + dvy*dvy + dvz*dvz);
if (dvcosmo < dvmin) {
dvmin=dvcosmo;
ivmin = i;
}
}
}
for (;;) {
double r2min = FLT_MAX;
q = NULL;
for (i=0;i<nSmooth;++i) {
if(TYPETest(nnList[i].p,TYPE_DELETED)) continue;
if(nnList[i].p->curlv()[2] == p->iOrder) continue;
/* don't choose a pre-accreted particle
but it can have been accreted by another BH */
#if 0
if(nnList[i].p->rung < s.iSinkCurrentRung) continue;
/* has to be on the right timestep */
#endif
r2 = nnList[i].fKey;
if(!s.bBHAccreteAll
&& r2 > 0.25*(nnList[i].p->fBall)*(nnList[i].p->fBall))
continue;
/* has to be nearby! */
rs = KERNEL(r2, nSmooth);
if(s.bBHMindv == 1)
weight = rs*pow(nnList[i].p->c()*nnList[i].p->c()
+(dvmin*dvmin),-1.5)/dCosmoDenFac;
else {
dvx = (-p->velocity[0]+nnList[i].p->vPred()[0])-aDot*nnList[i].dx.x;
dvy = (-p->velocity[1]+nnList[i].p->vPred()[1])-aDot*nnList[i].dx.y;
dvz = (-p->velocity[2]+nnList[i].p->vPred()[2])-aDot*nnList[i].dx.z;
dvcosmo = sqrt(dvx*dvx+dvy*dvy+dvz*dvz);
weight = rs*pow(nnList[i].p->c()*nnList[i].p->c()
+dvcosmo*dvcosmo,-1.5)/dCosmoDenFac; /* weight particles by mdot quantities */
/* cosmo factors put in 7/7/09 JMB */
}
if (weight > weat) {
r2min = r2; /* note r2min is not really the min r2 anymore */
weat = weight;
ieat = i;
q = nnList[i].p;
}
}
weat = -1e37; /* reset so other particles can be selected JMB 7/9/09 */
if(q == NULL) {
// dtEff = s.dSinkCurrentDelta; /**pow(0.5,iRung-smf->iSinkCurrentRung);*/
dtEff = RungToDt(dDelta, p->rung);
CkPrintf("BHSink %ld: WARNING!! Not enough edible particles. Time: %g Mass not eaten: %g \n",p->iOrder,dTime,mdotCurr*dtEff);
if(naccreted == 0) return;
else goto dofeedback;
}
else {
assert(q != NULL);
/* Timestep for accretion is larger of sink and victim timestep */
#if 0
iRung = q->rung;
if (iRung > p->rung) iRung = p->rung;
#endif
// dtEff = s.dSinkCurrentDelta;
dtEff = RungToDt(dDelta, p->rung);
dmq = mdotCurr*dtEff;
if (dmq < q->curlv()[0]) mdotCurr = 0.0; /* Original mass in q->curlv[0] */
else {
mdotCurr -= mdotCurr*(q->curlv()[0]/dmq); /* need an additional victim */
dmq = q->curlv()[0];
}
CkPrintf("BHSink %ld: Time: %g %ld dmq %g %g %g sharing \n",
p->iOrder,dTime,q->iOrder,dmq,q->curlv()[1],
p->dDeltaM());
ifMass = 1./(p->mass + dmq);
/* to record angular momentum JMB 11/9/10 */
CkPrintf("BHSink %ld: Gas: %ld dx: %g dy: %g dz: %g \n",
p->iOrder,q->iOrder,p->position[0]-q->position[0],
p->position[1]-q->position[1],
p->position[2]-q->position[2]);
CkPrintf("BHSink %ld: Gas: %ld dvx: %g dvy: %g dvz: %g \n",
p->iOrder,q->iOrder,p->velocity[0]-q->velocity[0],
p->velocity[1]-q->velocity[1],
p->velocity[2]-q->velocity[2]);
/* Adjust sink properties (conserving momentum etc...) */
p->position[0] = ifMass*(p->mass*p->position[0]
+dmq*q->position[0]);
p->position[1] = ifMass*(p->mass*p->position[1]
+dmq*q->position[1]);
p->position[2] = ifMass*(p->mass*p->position[2]
+dmq*q->position[2]);
p->velocity[0] = ifMass*(p->mass*p->velocity[0]
+dmq*q->velocity[0]);
p->velocity[1] = ifMass*(p->mass*p->velocity[1]
+dmq*q->velocity[1]);
p->velocity[2] = ifMass*(p->mass*p->velocity[2]
+dmq*q->velocity[2]);
p->treeAcceleration[0] = ifMass*(p->mass*p->treeAcceleration[0]
+dmq*q->treeAcceleration[0]);
p->treeAcceleration[1] = ifMass*(p->mass*p->treeAcceleration[1]
+dmq*q->treeAcceleration[1]);
p->treeAcceleration[2] = ifMass*(p->mass*p->treeAcceleration[2]
+dmq*q->treeAcceleration[2]);
p->fStarMetals() = ifMass*(p->mass*p->fStarMetals()+dmq*q->fMetals());
p->mass += dmq;
dm += dmq;
q->mass -= dmq;
naccreted += 1; /* running tally of how many are accreted JMB 10/23/08 */
CkPrintf("BHSink %ld: Time %g dist2 %ld %g gas smooth: %g eatenmass %g\n",
p->iOrder,dTime,q->iOrder,r2min,q->fBall*q->fBall,dmq);
if (q->mass <= dMinGasMass) {
if(!(TYPETest(q,TYPE_DELETED))) deleteParticle(q);
/* Particles are getting deleted twice, which is messing
up the bookkeeping. I think this will make them only be
deleted once. JMB 9/23/08*/
}
}
if (mdotCurr == 0.0) break;
}
}
/********************************************************************/
dofeedback:
if(s.dBHSinkFeedbackFactor != 0.0) { /* allows FB to be turned off JMB 7/20/09 */
dE = s.dBHSinkFeedbackFactor*dm; /* dE based on actual mass eaten */
// dtEff = s.dSinkCurrentDelta;
dtEff = RungToDt(dDelta, p->rung);
dmAvg = mdot*dtEff;
CkPrintf("BHSink %ld: Delta: %g Time: %.8f dm: %g dE %g\n",
p->iOrder,dtEff,dTime,dm,dE);
/* Recalculate Normalization */
ih2 = invH2(p);
f2h2 = p->fBall*p->fBall;
fNorm_new = 0.0;
fNorm = 0.5*M_1_PI*sqrt(ih2)*ih2;
for (i=0;i<nSmooth;++i) {
r2 = nnList[i].fKey*ih2;
rs = KERNEL(r2, nSmooth);
fNorm_new += rs;
rs *= fNorm;
}
/* spread feedback to nearest 32 particles across the
kernel. JMB 8/29/11 */
CkAssert(fNorm_new != 0.0);
fNorm_new = 1./fNorm_new;
counter = 0;
for (i=0;i<nSmooth;++i) {
r2 = nnList[i].fKey*ih2;
rs = KERNEL(r2, nSmooth);
fbweight = rs*fNorm_new;
nnList[i].p->u() += fbweight*dE;
nnList[i].p->uPred() += fbweight*dE;
counter ++;
/* now turn off cooling */
if( s.bBHTurnOffCooling)
nnList[i].p->fTimeCoolIsOffUntil()
= max(nnList[i].p->fTimeCoolIsOffUntil(),
dTime + RungToDt(dDelta, nnList[i].p->rung));
CkPrintf("BHSink %ld: Time %g FB Energy to %li dE %g tCoolOffUntil %g \n",
p->iOrder, dTime, nnList[i].p->iOrder,
fbweight*dE,nnList[i].p->fTimeCoolIsOffUntil());
if(s.bBHTurnOffCooling)
CkAssert(nnList[i].p->fTimeCoolIsOffUntil() > dTime);
nnList[i].p->dTimeFB() = dTime;
}
}
}
void BHAccreteSmoothParams::postTreeParticle(GravityParticle *p1)
{
/* Convert energy back to specific quantities (per mass)
because we are done with our conservative calculations */
if(TYPETest(p1, TYPE_GAS) && p1->mass != 0
&& !(TYPETest(p1,TYPE_DELETED))) {
p1->u() /= p1->mass;
p1->uPred() /= p1->mass;
}
}
int BHIdentifySmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
void BHIdentifySmoothParams::initSmoothCache(GravityParticle *p)
{
p->iEaterOrder() = 0;
}
void BHIdentifySmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if(p2->iEaterOrder > p1->iEaterOrder())
p1->iEaterOrder() = p2->iEaterOrder;
}
void BHIdentifySmoothParams::fcnSmooth(GravityParticle *p, int nSmooth,
pqSmoothNode *nnList)
{
/* Identify BHs for merging JMB */
GravityParticle *q = NULL;
int i;
double deltaa, deltar, deltav;
double aFac = a;
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
if(nnList[i].fKey > 4.0*p->soft*p->soft || q->iOrder == p->iOrder)
continue;
/* are they close together? (within TWO softenings JMB 10/21/10*/
/* don't include yourself! JMB 12/11/08 */
deltaa=sqrt(sqr(p->treeAcceleration[0] - q->treeAcceleration[0])
+ sqr(p->treeAcceleration[1] - q->treeAcceleration[1])
+ sqr(p->treeAcceleration[2] - q->treeAcceleration[2]));
deltar=sqrt(sqr(p->position[0] - q->position[0])
+ sqr(p->position[1] - q->position[1])
+ sqr(p->position[2] - q->position[2]));
deltav=sqrt(sqr(p->velocity[0] - q->velocity[0])
+ sqr(p->velocity[1] - q->velocity[1])
+ sqr(p->velocity[2] - q->velocity[2]));
if ( deltaa*deltar*aFac < 0.5*deltav*deltav ) continue;
/* Selects other BH particles that are
* within the criteria
* delta_a*delta_r < .5*delta_v^2
AND
* within the softening */
if(p->iOrder < q->iOrder) {
if(p->iOrder < q->iEaterOrder() || q->iEaterOrder() == 0) {
q->iEaterOrder() = p->iOrder;
CkPrintf("BHSink MergeID %ld will be eaten by %ld \n",
q->iOrder,p->iOrder);
/* iEaterOrder is the place holder for the iord of the
eating black hole. the victim remembers who eats
it for the merging step which is next. JMB 4/30/09 */
}
}
}
}
int BHSinkMergeSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
void BHSinkMergeSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
/* Makes BHs merge together, based on the criteria
set in BHSinkIdentify above */
GravityParticle *q = NULL;
double ifMass;
int i;
/* for kicks */
const double A = 12000.0;
const double B = -0.93;
const double H = 7300.0;
const double K = 60000.0; /* BH kick parameters in km/s */
double vkick, vm, vpar, vkx, vky, vkz;
double vperp, mratio, mfactor, cosine, spin1, spin2, angle1, angle2, xi;
double xrand, yrand, zrand;
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
/* Give mass to higher iOrd particle, delete lower. */
if(q->iEaterOrder() == p->iOrder) {
/* q has already been identified as a victim of p */
ifMass = 1./(p->mass + q->mass);
mratio = p->mass/q->mass;
if(mratio > 1.0) mratio = q->mass/p->mass;
/* Adjust sink properties (conserving momentum etc...) */
p->position[0] = ifMass*(p->mass*p->position[0]
+q->mass*q->position[0]);
p->position[1] = ifMass*(p->mass*p->position[1]
+q->mass*q->position[1]);
p->position[2] = ifMass*(p->mass*p->position[2]
+q->mass*q->position[2]);
p->velocity[0] = ifMass*(p->mass*p->velocity[0]
+q->mass*q->velocity[0]);
p->velocity[1] = ifMass*(p->mass*p->velocity[1]
+q->mass*q->velocity[1]);
p->velocity[2] = ifMass*(p->mass*p->velocity[2]
+q->mass*q->velocity[2]);
p->treeAcceleration[0] = ifMass*(p->mass*p->treeAcceleration[0]
+q->mass*q->treeAcceleration[0]);
p->treeAcceleration[1] = ifMass*(p->mass*p->treeAcceleration[1]
+q->mass*q->treeAcceleration[1]);
p->treeAcceleration[2] = ifMass*(p->mass*p->treeAcceleration[2]
+q->mass*q->treeAcceleration[2]);
p->mass += q->mass;
/**** Gravitational Recoil ****/
if(s.bDoBHKick == 1
&& (fabs(-1.0*p->fTimeForm() - dTime) > s.dDeltaStarForm)) {
/* Turn off recoil if BH has just formed. Ejecting them
immediately is not helpful. JMB 8/5/09 */
mfactor = pow(mratio,2)/pow((1+mratio),5);
vm = A * (1.0-mratio) * mfactor * (1.0 + B*mratio/pow((1.0+mratio),2));
cosine = (rand()/((double) RAND_MAX )) * 2.*M_PI; /* angle */
spin1 = (rand()/((double) RAND_MAX )); /* spins of BHs 1 and 2 */
spin2 = (rand()/((double) RAND_MAX ));
angle1 = (rand()/((double) RAND_MAX )) * 2.*M_PI; /* orientations of BH spins */
angle2 = (rand()/((double) RAND_MAX )) * 2.*M_PI;
xi = 1.5358897; /* 88.0 * 2. * 3.14159 / 360.0 */
/* in Campanelli et al xi = 88 degrees */
vperp = H * mfactor * (spin2*cos(angle2) - mratio*spin1*cos(angle1)); /* check sin cos! */
vpar = K * cos(cosine) * mfactor * (spin2*sin(angle2) - mratio*spin1*sin(angle1));
vkick = sqrt ( (vm + vperp*cos(xi))*(vm + vperp*cos(xi)) + (vperp*sin(xi))*(vperp*sin(xi)) + vpar*vpar);
vkick = vkick / s.dKmPerSecUnit;
/* comoving coords are important JMB 5/15/09 */
vkick = vkick*a;
/* random direction */
for(;;){
xrand = (rand()/((double) RAND_MAX )) * 2.0 - 1.0;
yrand = (rand()/((double) RAND_MAX )) * 2.0 - 1.0;
zrand = (rand()/((double) RAND_MAX )) * 2.0 - 1.0;
if(sqrt(xrand*xrand+yrand*yrand+zrand*zrand) < 1.0) break;
}
vkx = vkick * xrand / sqrt(xrand*xrand+yrand*yrand+zrand*zrand);
vky = vkick * yrand / sqrt(xrand*xrand+yrand*yrand+zrand*zrand);
vkz = vkick * zrand / sqrt(xrand*xrand+yrand*yrand+zrand*zrand);
p->velocity[0] += vkx;
p->velocity[1] += vky;
p->velocity[2] += vkz;
}
else vkick = 0.0;
CkPrintf("BHSink Merge %ld eating %ld Time %g kick velocity %g mass ratio %g \n",
p->iOrder,q->iOrder,dTime,vkick,mratio);
if(q->iEaterOrder() > 0 && !(TYPETest(q,TYPE_DELETED))) {
CkPrintf("BHSink Merge %ld eaten Time %g \n", q->iOrder,dTime);
deleteParticle(q);
}
}
}
}
void BHSinkMergeSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if (!(TYPETest( p1, TYPE_DELETED )) && TYPETest( p2, TYPE_DELETED ) ) {
deleteParticle(p1);
}
}
int SinkAccreteTestSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
static inline double & fBindingEnergy(GravityParticle *a) {return a->curlv()[0];}
static inline double & fBindingEnergy(ExternalSmoothParticle *a)
{return a->curlv[0];}
/* Indicator for r,v,a update */
#define bRVAUpdate(_a) ((_a)->curlv()[2])
void SinkAccreteTestSmoothParams::initSmoothCache(GravityParticle *p)
{
fBindingEnergy(p) = FLT_MAX;
iOrderSink(p) = -1;
#ifdef SINKING
if (TYPETest(p, TYPE_SINKING)) {
iOrderSink(p) = p->iSinkingOnto;
}
#endif
}
void SinkAccreteTestSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
/* Particle p1 belongs to sink iOrderSink(p1) initially but will
switch to iOrderSink(p2) if more bound to that sink */
#ifdef SINKING
if (TYPETest(p1, TYPE_SINKING)) return;
#endif
if (fBindingEnergy(p2) < fBindingEnergy(p1)) {
fBindingEnergy(p1) = fBindingEnergy(p2);
iOrderSink(p1) = iOrderSink_Ex(p2);
}
if (TYPETest( p2, TYPE_NEWSINKING)) TYPESet( p1, TYPE_NEWSINKING);
}
#ifdef SINKING
#define TRUEMASS(p__) (p__)->fTrueMass
#else
#define TRUEMASS(p__) (p__)->mass
#endif
void SinkAccreteTestSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
int i;
double dSinkRadius2 = s.dSinkRadius*s.dSinkRadius,
EBO,Eq,r2,dvx,dv2;
GravityParticle *q;
#ifdef SINKING
assert(TYPETest(p, TYPE_SINK));
assert(!TYPETest(p, TYPE_SINKING));
#endif
/* G = 1
p is sink particle
q is gas particle */
if (s.dSinkBoundOrbitRadius > 0)
EBO = -0.5*TRUEMASS(p)/s.dSinkBoundOrbitRadius;
else
EBO = FLT_MAX;
for (i=0;i<nSmooth;++i) {
r2 = nnList[i].fKey;
if (r2 > 0 && r2 <= dSinkRadius2) {
q = nnList[i].p;
#ifdef SINKING
if (TYPETest(q, TYPE_SINKING)) {
FLOAT r0 = q->rSinking0Mag;
FLOAT r2 = r0 + q->vSinkingr0*(smf->dTime-q->fSinkingTime);
continue;
}
#endif
if (TYPETest( q, TYPE_GAS )
&& (q->rung >= s.iSinkCurrentRung
|| TYPETest(q, TYPE_NEWSINKING))) {
dvx = p->velocity[0]-q->velocity[0];
dv2 = dvx*dvx;
dvx = p->velocity[1]-q->velocity[1];
dv2 += dvx*dvx;
dvx = p->velocity[2]-q->velocity[2];
dv2 += dvx*dvx;
Eq = -TRUEMASS(p)/sqrt(r2) + 0.5*dv2;
if (s.bSinkThermal) Eq+= q->u();
/* Being labelled NEWSINKING forces accretion to somebody even if unbound */
if (Eq < EBO || r2 < sqr(s.dSinkMustAccreteRadius)
|| TYPETest(q, TYPE_NEWSINKING)) {
if (Eq < fBindingEnergy(q)) {
fBindingEnergy(q) = Eq;
iOrderSink(q) = p->iOrder; /* Particle q belongs to sink p */
TYPESet(q, TYPE_NEWSINKING);
}
}
}
}
}
}
int SinkingAverageSmoothParams::isSmoothActive(GravityParticle *p)
{
return (TYPETest(p, TYPE_NEWSINKING));
}
void SinkingAverageSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
#ifdef SINKING
int i,j;
PARTICLE *q;
double dv[3],wt,ih2,rs,r2,norm;
assert(TYPETest(p,TYPE_NEWSINKING));
wt = 0;
for (j=0;j<3;j++) dv[j]=0;
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
ih2 = 4.0/BALL2(p);
r2 = nnList[i].fDist2*ih2;
KERNEL(rs,r2);
rs*= q->fMass;
wt+=rs;
for (j=0;j<3;j++) dv[j]+=rs*(q->v[j]-p->v[j]);
}
norm=1./wt;
/* smoothed velocity */
for (j=0;j<3;j++) {
p->vSinkingTang0Unit[j] = dv[j]*norm;
}
#endif
}
int SinkAccreteSmoothParams::isSmoothActive(GravityParticle *p)
{
if(p->rung < activeRung)
return 0;
return (TYPETest(p, TYPE_SINK));
}
void SinkAccreteSmoothParams::initSmoothCache(GravityParticle *p)
{ /* cached copies only */
#ifdef SINKING
if (TYPETest( ((PARTICLE *) p), TYPE_SINKING )) {
bRVAUpdate(p) = 0; /* Indicator for r,v,a update */
}
#endif
}
void SinkAccreteSmoothParams::combSmoothCache(GravityParticle *p1,
ExternalSmoothParticle *p2)
{
if (!(TYPETest( p1, TYPE_DELETED )) &&
TYPETest( p2, TYPE_DELETED ) ) {
p1->mass = p2->mass;
TYPEReset(p1,TYPE_SINKING);
deleteParticle(p1);
}
#ifdef SINKING
else if (TYPETest(p2, TYPE_SINKING ) && !TYPETest( p1, TYPE_SINKING )) {
GravityParticle *p=p1,*q=p2;
int j;
#ifdef SINKDBG
printf("SINKING Sinking particle on other processor -- sinking properties combined %ld %g\n",p->iOrder,p->fMetals);
#endif
bRVAUpdate(p) = 1;
assert(!TYPETest(p,TYPE_SINK));
assert(!TYPETest(q,TYPE_SINK));
p->fTrueMass = q->fTrueMass;
TYPESet(p, TYPE_SINKING);
p->fMetals = -1; /* HACK -- flag in output for sinking state */
p->iSinkingOnto = q->iSinkingOnto;
p->dt = q->dt;
p->rung = q->iRung;
p->fSinkingTime = q->fSinkingTime;
p->vSinkingr0 = q->vSinkingr0;
p->rSinking0Mag = q->rSinking0Mag;
p->vSinkingTang0Mag = q->vSinkingTang0Mag;
assert(q->vSinkingr0 < 0);
for (j=0;j<3;j++) {
p->rSinking0Unit[j] = q->rSinking0Unit[j];
p->vSinkingTang0Unit[j] = q->vSinkingTang0Unit[j];
}
}
if (TYPETest( p1, TYPE_SINKING )) {
#ifdef SINKINGAVERAGE
TYPEReset(p1,TYPE_NEWSINKING);
#endif
if (bRVAUpdate(p2)) {
GravityParticle *p=p1,*q=p2;
p->r[0] = q->r[0];
p->r[1] = q->r[1];
p->r[2] = q->r[2];
p->v[0] = q->v[0];
p->v[1] = q->v[1];
p->v[2] = q->v[2];
p->vPred[0] = q->vPred[0];
p->vPred[1] = q->vPred[1];
p->vPred[2] = q->vPred[2];
p->a[0] = q->a[0];
p->a[1] = q->a[1];
p->a[2] = q->a[2];
#ifdef SINKDBG
if (((PARTICLE *)p1)->iOrder == 80) printf("FORCESINKACCRETECOMB %d with %d \n",-1,((PARTICLE *)p1)->iOrder);
#endif
}
}
#endif
}
/*#define DBGIORDER 2000000 */
#define DBGIORDER -1
void SinkAccreteSmoothParams::fcnSmooth(GravityParticle *p,int nSmooth,
pqSmoothNode *nnList)
{
int i,j,bEat=0;
double ifMass;
GravityParticle *q;
#ifdef SINKING
double drSink[3];
double r2;
for (j=0;j<3;j++) drSink[j] = -p->position[j];
#endif
#ifdef STARSINK
/* Convert to total angular momentum for additions */
SINK_Lx(p) += TRUEMASS(p)*(p->position[1]*p->velocity[2]
- p->position[2]*p->velocity[1]);
SINK_Ly(p) += TRUEMASS(p)*(p->position[2]*p->velocity[0]
- p->position[0]*p->velocity[2]);
SINK_Lz(p) += TRUEMASS(p)*(p->position[0]*p->velocity[1]
- p->position[1]*p->velocity[0]);
#endif
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
if ( iOrderSink(q) == p->iOrder) {
#ifndef SINKING
double mp,mq,rx,ry,rz,vx,vy,vz;
mp = p->mass; mq = q->mass;
rx = q->position[0]; ry = q->position[1]; rz = q->position[2];
vx = q->velocity[0]; vy = q->velocity[1]; vz = q->velocity[2];
ifMass = 1./(mp+mq);
for (j=0;j<3;j++) {
p->position[j] = ifMass*(mp*p->position[j]
+mq*q->position[j]);
p->velocity[j] = ifMass*(mp*p->velocity[j]
+mq*q->velocity[j]);
p->treeAcceleration[j] = ifMass*(mp*p->treeAcceleration[j]
+mq*q->treeAcceleration[j]);
}
#ifdef STARSINK
SINK_Lx(p) += mq*(ry*vz - rz*vy); /* add to total L */
SINK_Ly(p) += mq*(rz*vx - rx*vz);
SINK_Lz(p) += mq*(rx*vy - ry*vx);
#endif
/* Add internal energy to sink particle if it is gas*/
if (p->isGas()) p->u() += ifMass*(mp*p->u()+mq*q->u());
p->mass += mq;
bEat = 1;
assert(q->mass != 0);
q->mass = 0;
deleteParticle(q);
#endif
#ifdef SINKING
assert(TYPETest(q, TYPE_GAS));
assert(!TYPETest(q, TYPE_SINK));
r2 = nnList[i].fKey;
if (r2 < sinks.dSinkMustAccreteRadius*sinks.dSinkMustAccreteRadius ||
(TYPETest(q, TYPE_SINKING) &&
q->rSinking0Mag + q->vSinkingr0*(dTime-q->fSinkingTime) < smf->dSinkMustAccreteRadius)) {
#ifdef SINKDBG
if (q->iOrder == 55) printf("FORCESINKACCRETE0 %ld with %ld \n",p->iOrder,q->iOrder);
#endif
if (!TYPETest(q, TYPE_SINKING)) {
double mp,mq,rx,ry,rz,vx,vy,vz;
mp = p->fTrueMass; mq = q->fMass;
rx = q->r[0]; ry = q->r[1]; rz = q->r[2];
vx = q->v[0]; vy = q->v[1]; vz = q->v[2];
/* Do a standard accrete -- no sinking stage */
ifMass = 1./(mp+mq);
for (j=0;j<3;j++) {
p->r[j] = ifMass*(mp*p->r[j]+mq*q->r[j]);
p->v[j] = ifMass*(mp*p->v[j]+mq*q->v[j]);
p->a[j] = ifMass*(mp*p->a[j]+mq*q->a[j]);
}
#ifdef STARSINK
SINK_Lx(p) += mq*(ry*vz - rz*vy); /* add to total L */
SINK_Ly(p) += mq*(rz*vx - rx*vz);
SINK_Lz(p) += mq*(rx*vy - ry*vx);
#endif
p->u += ifMass*(mp*p->u+mq*q->u);
p->fMass += mq;
p->fTrueMass += q->fMass;
if (p->iOrder==DBGIORDER) printf("SINKING %ld Direct Delete of %ld +%g %g %g\n",p->iOrder,q->iOrder,q->fMass,p->fMass,p->fTrueMass);
bEat = 1;
assert(q->fMass != 0);
q->fMass = 0;
pkdDeleteParticle(smf->pkd, q);
}
else if (q->iSinkingOnto == p->iOrder) {
/* Already sinking -- now do final accrete */
p->fMass += q->fMass; /* Note: p->fTrueMass NOT modified */
p->iSinkingOnto--; /* One less sinking onto this sink */
if (p->iOrder==DBGIORDER) printf("SINKING %ld Final Delete of %ld +%g %g %g %d\n",p->iOrder,q->iOrder,q->fMass,p->fMass,p->fTrueMass,p->iSinkingOnto);
bEat = 1;
assert(q->fMass != 0);
q->fMass = 0;
TYPEReset(q,TYPE_SINKING); /* Don't want it to participate later */
pkdDeleteParticle(smf->pkd, q);
}
}
else if (!TYPETest(q, TYPE_SINKING)) {
FLOAT mp,mq,rx,ry,rz,vx,vy,vz;
mp = p->fTrueMass; mq = q->fMass;
rx = q->r[0]; ry = q->r[1]; rz = q->r[2];
vx = q->v[0]; vy = q->v[1]; vz = q->v[2];
/* Enter sinking stage */
assert(r2 < smf->dSinkRadius*smf->dSinkRadius);
/* If this code is active, ignore non-infalling particles */
/* vr = 0;
for (j=0;j<3;j++) {
vr += (q->r[j]-p->r[j])*(q->v[j]-p->v[j]);
}
if (vr >= 0) continue;*/
/* All force, velocity, position info associated
with particle now belong to sink instead.
Note: If we want accurate estimate of L for the sink
we should adjust the sink L now */
/* Everything short of eating the particle
If many additional particles are added sink position could vary wildly
-- defer trajectory calculation until after */
ifMass = 1./(p->fTrueMass+q->fMass);
for (j=0;j<3;j++) {
p->r[j] = ifMass*(mp*p->r[j]+q->fMass*q->r[j]);
p->v[j] = ifMass*(mp*p->v[j]+q->fMass*q->v[j]);
p->a[j] = ifMass*(mp*p->a[j]+q->fMass*q->a[j]);
}
#ifdef STARSINK
SINK_Lx(p) += mq*(ry*vz - rz*vy); /* Add to total L */
SINK_Ly(p) += mq*(rz*vx - rx*vz);
SINK_Lz(p) += mq*(rx*vy - ry*vx);
#endif
p->u += ifMass*(mp*p->u+mq*q->u);
p->fTrueMass += q->fMass;
p->iSinkingOnto++; /* One more sinking onto this sink */
q->fTrueMass = 0; /* Sinking particles retain mass to act as ghost particles for forces g+SPH */
assert(!TYPETest(q,TYPE_SINK));
assert(TYPETest(q,TYPE_NEWSINKING));
TYPESet(q, TYPE_SINKING);
q->fMetals = -1; /* HACK -- flag for sinking state */
q->iSinkingOnto = p->iOrder;
q->dt = p->dt;
q->iRung = p->iRung;
q->fSinkingTime = smf->dTime;
q->vSinkingr0 = 1e30;
if (p->iOrder==DBGIORDER) printf("SINKING %ld Add Sinking %ld +%g %g %g %d\n",p->iOrder,q->iOrder,q->fMass,p->fMass,p->fTrueMass,p->iSinkingOnto);
bEat = 1;
}
#endif /*SINKING*/
}
}
#ifdef STARSINK
/* Store Internal angular momentum only as it is invariant under motions */
SINK_Lx(p) -= TRUEMASS(p)*(p->position[1]*p->velocity[2]
- p->position[2]*p->velocity[1]);
SINK_Ly(p) -= TRUEMASS(p)*(p->position[2]*p->velocity[0]
- p->position[0]*p->velocity[2]);
SINK_Lz(p) -= TRUEMASS(p)*(p->position[0]*p->velocity[1]
- p->position[1]*p->velocity[0]);
#endif
#ifdef SINKING
if (bEat) {
for (j=0;j<3;j++) drSink[j] += p->r[j]; /* If sink ate it moved: old sinking particles must follow */
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
if (TYPETest( q, TYPE_SINKING ) && q->iSinkingOnto == p->iOrder) {
if (TYPETest(q,TYPE_NEWSINKING)) {
double r0, dx[3], dv[3], vr, vrraw, vrmin, vt, norm, dot;
/* Initialize sinking trajectory */
r0 = 0;
for (j=0;j<3;j++) {
dx[j] = (q->r[j]-p->r[j]);
r0 += dx[j]*dx[j];
}
r0 = sqrt(r0);
norm = 1/r0;
for (j=0;j<3;j++) dx[j] *= norm;
vr = 0;
for (j=0;j<3;j++) {
dv[j] = q->v[j]-p->v[j]; /* dv = raw relative velocity */
vr += dv[j]*dx[j];
} /* raw infall v */
#ifdef SINKINGAVERAGE
vrraw = vr;
/* See if smoothed velocity is still infalling
q->vSinkingTang0Unit temporarily contains mean of (q'-q) */
vr = 0;
for (j=0;j<3;j++) {
dv[j] = q->vSinkingTang0Unit[j]+q->v[j]-p->v[j]; /* dv now smoothed relative vel */
vr += dv[j]*dx[j];
}
if (p->iOrder==DBGIORDER) printf("SINKING %ld Add Sinking Trajectory vr %ld %g %g %g\n",p->iOrder,q->iOrder,vr,vrraw,-q->c);
/* printf("Infall vr avg %g ",vr);*/
#endif
vt = 0;
for (j=0;j<3;j++) {
dv[j] -= vr*dx[j]; /*subtract off radial part*/
vt += dv[j]*dv[j];
}
vt = sqrt(vt); /* Averaged tangential motion */
#ifdef SINKINGAVERAGE
if (vr >= 0) vr=vrraw;
#endif
vrmin = -q->c; /* Bondi-Hoyle minimum inflow speed -- perhaps not ideal for rotating case */
if (vr > vrmin) vr = vrmin;
/* printf("min %g : %g vt %g \n",-vrmin,vr,vt);*/
q->rSinking0Mag = r0;
if (r0 > smf->dSinkRadius) printf("WARNING: %li outside sink r=%g %g %g\n",q->iOrder,r0,sqrt(nnList[i].fDist2),smf->dSinkRadius);
q->vSinkingr0 = vr;
q->vSinkingTang0Mag = vt;
norm = 1/(vt+1e-37);
dot = 0;
for (j=0;j<3;j++) {
q->rSinking0Unit[j] = dx[j];
q->vSinkingTang0Unit[j] = dv[j]*norm;
dot += q->rSinking0Unit[j]*q->vSinkingTang0Unit[j]; /* sanity check */
q->a[j] = p->a[j];
q->v[j] = p->v[j]; /* move with sink now */
}
assert(fabs(dot) < 1e-4);
}
else {
assert(q->vSinkingr0 < 0);
for (j=0;j<3;j++) {
q->r[j] += drSink[j];
q->a[j] = p->a[j];
q->v[j] = p->v[j];
}
}
}
}
}
/* Make sure sinking particles are still on their trajectories */
for (i=0;i<nSmooth;++i) {
q = nnList[i].p;
if (TYPETest(q,TYPE_SINKING) && q->iSinkingOnto == p->iOrder) {
FLOAT r0 = q->rSinking0Mag;
FLOAT r2 = r0 + q->vSinkingr0*(smf->dTime-q->fSinkingTime);
FLOAT thfac, th2, costh2, sinth2, dr2, sqr02;
FLOAT dr[3];
assert(q->vSinkingr0 < 0);
if (r2 < 0.1*r0) r2 = 0.1*r0; /* HACK */
thfac = q->vSinkingTang0Mag*2/(q->vSinkingr0);
sqr02 = sqrt(r0/r2);
th2 = thfac*(1-sqr02);
costh2 = cos(th2);
sinth2 = sin(th2);
for (j=0;j<3;j++) {
q->r[j] = p->r[j]+r2*costh2*q->rSinking0Unit[j]+r2*sinth2*q->vSinkingTang0Unit[j];
q->vPred[j] = p->v[j]+q->vSinkingTang0Mag*sqr02*
(-sinth2*q->rSinking0Unit[j]+costh2*q->vSinkingTang0Unit[j])
+q->vSinkingr0*(costh2*q->rSinking0Unit[j]+sinth2*q->vSinkingTang0Unit[j]);
}
bRVAUpdate(q) = 1; /* Indicator for r,v,a update */
#ifdef SINKDBG
if (q->iOrder == 160460) printf("SINKINGFORCESHARE %ld %g %g %g %g %g %g %g %g %g (%g)\n",q->iOrder,smf->dTime,r2,sqrt(pow(q->r[0]-p->r[0],2.)+pow(q->r[1]-p->r[1],2.)+pow(q->r[2]-p->r[2],2.)), q->r[0]-p->r[0],q->r[1]-p->r[1],q->r[2]-p->r[2],q->vPred[0]-p->v[0],q->vPred[1]-p->v[1],q->vPred[2]-p->v[2],(q->vPred[0]-p->v[0])*(q->r[0]-p->r[0])+(q->vPred[1]-p->v[1])*(q->r[1]-p->r[1])+(q->vPred[2]-p->v[2])*(q->r[2]-p->r[2]));
#endif
TYPEReset(q,TYPE_NEWSINKING);
}
}
#endif
}
