//SIDM- self interacting dark matter by Alexander Fry
#include <math.h>
#ifdef HAVE_VALUES_H
#include <values.h>
#else
#include <float.h>
#endif
#include "ParallelGravity.h"
#include "smooth.h"
#include "SIDM.h"
///
/// @brief Main method to perform Self Interacting Dark Matter interactions
/// @param dTime current simulation time
/// @param dDelta timestep over which to calculate interactions
/// @param activeRung timestep rung corresponding to dDelta
///
void Main::doSIDM(double dTime,double dDelta, int activeRung) {
if (param.iSIDMSelect!=0) {
double stime = CkWallTimer();
if(verbosity > 0) CkPrintf("SIDM interactions ... ");
SIDMSmoothParams pSIDM(TYPE_DARK, activeRung, param.csm, dTime,param.dSIDMSigma, param.dSIDMVariable,param.iSIDMSelect,param.dDelta );
treeProxy.startSmooth(&pSIDM, 1, param.nSmooth, 0, CkCallbackResumeThread());
if(verbosity > 0) CkPrintf("took %g seconds\n", CkWallTimer() - stime);
}
}
int SIDMSmoothParams::isSmoothActive(GravityParticle *p) {
return (TYPETest(p, iType));
}
void SIDMSmoothParams::initSmoothParticle(GravityParticle *p) {
}
void SIDMSmoothParams::initTreeParticle(GravityParticle *p) {
}
void SIDMSmoothParams::initSmoothCache(GravityParticle *p) {
p->treeAcceleration = p->velocity;
}
void SIDMSmoothParams::combSmoothCache(GravityParticle *p1, ExternalSmoothParticle *p2) {
Vector3D<double> deltav; //vector template notation
deltav=p2->velocity - p2->treeAcceleration;
p1->velocity += deltav;
}
void SIDMSmoothParams::fcnSmooth(GravityParticle *p, int nSmooth, pqSmoothNode *nnList) {
GravityParticle *q;
double fNorm,ih2,r2;
double ran, probability,aDot, dvx,dvy,dvz,dvcosmo;
int i;
double mu, pxcm,pycm,pzcm,pcm0,ux,uy,uz;
double norm,vxcm,vycm,vzcm, m1,m2;
double uvar,vvar;
double sigma_classic, beta, term; //velocity depedence terms
double Sigma = 0.0;
aDot=a*H;
ih2 = invH2(p); //smoothing length, 1/h^2
fNorm = M_1_PI*ih2*sqrt(ih2);
for (i=0;i<nSmooth;++i)
{
q = nnList[i].p;
if (q->iOrder != p->iOrder) //don't interact with self
{
ran=rand()/((double)RAND_MAX); //random number on [0,1]
dvx = (-p->velocity.x + q->velocity.x)/a - aDot*nnList[i].dx.x; //v or vpred?
dvy = (-p->velocity.y + q->velocity.y)/a - aDot*nnList[i].dx.y;
dvz = (-p->velocity.z + q->velocity.z)/a - aDot*nnList[i].dx.z;
dvcosmo = sqrt(dvx*dvx + dvy*dvy + dvz*dvz); //relative velocity between particles
r2 = nnList[i].fKey*ih2;
if (iSIDMSelect==1)
{
Sigma=dSIDMSigma;
}
else if (iSIDMSelect==2) //classical cross section dSIDMVariable is break in km/s
{
beta=M_PI*dSIDMVariable*dSIDMVariable/(dvcosmo*dvcosmo);
if (beta < .1){
sigma_classic=(4.0*M_PI/ 22.7 ) * beta *beta* log(1.0+(1.0/beta));
}
if ((beta >=.1 ) && (beta<1000)) {
sigma_classic=(8.0*M_PI/ 22.7 ) * beta *beta /( 1.0+1.5* pow(beta,1.65) );
}
if (beta >= 1000.0 ){
term= log(beta)+ 1.0 - .5*(1.0/log(beta));
sigma_classic=(M_PI/ 22.7 ) *term *term;
}
Sigma=sigma_classic*dSIDMSigma; //convert to simulation units
}
else if (iSIDMSelect==3) //resonant cross section, dSIDMVariable is exponent value
{
Sigma=dSIDMSigma*pow(dvcosmo,dSIDMVariable);
}
else
CkAbort("SIDM done with unknown cross section type (iSIDMSelect");
//density in physical units? fNorm*KERNEL( r2 )*(p->mass)/(a*a*a)
probability=Sigma*dvcosmo*dDelta*fNorm*KERNEL( r2, nSmooth )*(q->mass)/(a*a*a*2.0);
//if ( ran>0.999999) {
// CkPrintf("SIDM Diagnostics: %g \n",probability);
// CkPrintf("iOrder: %i, dDelta: %g, dvcosmo: %g, dSIDMSigma: %g \n",p->iOrder,dDelta,dvcosmo, dSIDMSigma);
// CkPrintf("fnorm: %g, kern: %g, r2: %g, q->mass: %g, aaa: %g \n",fNorm,KERNEL( r2 ),r2,(q->mass),(a*a*a*2.0));
// CkPrintf("probability: %g, dvcosmo: %g, sigma: %g \n", probability,dvcosmo,Sigma);
// CkPrintf("iSIDMSelect: %i \n",(int)iSIDMSelect);
// }
if (probability > .1 && ran>0.99999) {
CkPrintf("SIDM Warning! The probability is rather large: %g \n",probability);
CkPrintf("iOrder: %ld, dDelta: %g, dvcosmo: %g, Sigma: %g \n",p->iOrder,dDelta,dvcosmo, Sigma);
CkPrintf("fnorm: %g, kern: %g, r2: %g, q->mass: %g, aaa: %g \n",fNorm,KERNEL( r2, nSmooth ),r2,(q->mass),(a*a*a*2.0));
}
if (probability > ran) {
#ifdef SIDMINTERACT
p->iNSIDMInteractions += 1;
#endif
m1=p->mass;
m2=q->mass;
mu=m1*m2/(m1+m2);
pxcm = -mu*dvx; //momentum COM frame where vx1 is 0
pycm = -mu*dvy;
pzcm = -mu*dvz;
vxcm=-pxcm/m1; //velocity in COM frame where vx1 0
vycm=-pycm/m1;
vzcm=-pzcm/m1;
pcm0=sqrt(pxcm*pxcm+pycm*pycm+pzcm*pzcm);
norm=666;
while (norm>1.0){ //unit sphere point picking (Marsaglia 1972)
uvar=2.0*(rand()/(double)RAND_MAX)-1.0; //#random number on [-1,1]
vvar=2.0*(rand()/(double)RAND_MAX)-1.0;
norm=(uvar*uvar+vvar*vvar);
}
ux=2.0*uvar*sqrt(1.0-uvar*uvar-vvar*vvar);
uy=2.0*vvar*sqrt(1.0-uvar*uvar-vvar*vvar);
uz=1.0-2.0*(uvar*uvar+vvar*vvar);
pxcm=pcm0*ux;
pycm=pcm0*uy;
pzcm=pcm0*uz;
p->velocity.x+=a*(pxcm/m1);
p->velocity.y+=a*(pycm/m1);
p->velocity.z+=a*(pzcm/m1);
q->velocity.x+=a*(-pxcm/m2);
q->velocity.y+=a*(-pycm/m2);
q->velocity.z+=a*(-pzcm/m2);
//to transform back to the lab frame we must add back in the COM velocity
p->velocity.x+=a*(vxcm);
p->velocity.y+=a*(vycm);
p->velocity.z+=a*(vzcm);
q->velocity.x+=a*(vxcm-dvx);
q->velocity.y+=a*(vycm-dvy);
q->velocity.z+=a*(vzcm-dvz);
}
}
}
}