Homec4science
Diffusion Gear (master) Browse

Gear/src/sticky.c
73a76d71691cmaster

Last Change
Raw File

sticky.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 PERIODIC
static double boxSize, boxHalf;
#ifdef LONG_X
static double boxSize_X, boxHalf_X;
#else
#define boxSize_X boxSize
#define boxHalf_X boxHalf
#endif
#ifdef LONG_Y
static double boxSize_Y, boxHalf_Y;
#else
#define boxSize_Y boxSize
#define boxHalf_Y boxHalf
#endif
#ifdef LONG_Z
static double boxSize_Z, boxHalf_Z;
#else
#define boxSize_Z boxSize
#define boxHalf_Z boxHalf
#endif
#endif
#ifdef MULTIPHASE
static int Ncz;
static int Ncxy;
static int Ncxyz;
static double hubble_a, atime, hubble_a2, hubble_a_a, ainv, fac_mu, fac_vsic_fix, a3inv, fac_egy;
/*! \file sticky.c
* \brief Compute sticky collisions
*
*/
/*! init sticky
*
*/
void init_sticky()
{
Ncz = All.StickyGridNz;
Ncxy = All.StickyGridNx*All.StickyGridNy;
Ncxyz = All.StickyGridNx*All.StickyGridNy*All.StickyGridNz;
}
/*! This function is the driver routine for the calculation of sticky collisions
*/
void sticky()
{
double t0, t1;
double dt;
int Tis,Tic;
int do_sticky=0;
t0 = second(); /* measure the time for the full chimie computation */
do_sticky=0;
if (All.ComovingIntegrationOn)
{
Tis = log(All.StickyLastCollisionTime/All.TimeBegin) / All.Timebase_interval;
Tic = All.Ti_Current;
dt = get_cosmictime_difference(Tis,Tic);
if(dt >= All.StickyCollisionTime)
do_sticky=1;
}
else
{
if((All.Time - All.StickyLastCollisionTime) >= All.StickyCollisionTime)
do_sticky=1;
}
if(do_sticky)
{
sticky_compute_energy_kin(1);
//sticky_collisions(); /* old implementation */
if(ThisTask==0)
printf("sticky first loop\n");
sticky_collisions2(1);
if(ThisTask==0)
printf("sticky second loop\n");
sticky_collisions2(2);
sticky_compute_energy_kin(2);
//All.StickyLastCollisionTime += All.StickyCollisionTime;
All.StickyLastCollisionTime = All.Time;
}
t1 = second();
All.CPU_Sticky += timediff(t0, t1);
}
/*! This function is used as a comparison kernel in a sort routine.
*/
int compare_sticky_index(const void *a, const void *b)
{
if(((struct Sticky_index *) a)->CellIndex < (((struct Sticky_index *) b)->CellIndex))
return -1;
if(((struct Sticky_index *) a)->CellIndex > (((struct Sticky_index *) b)->CellIndex))
return +1;
return 0;
}
/*! This function returns the first particle that may sticky-interact with
* the given particle. This routine use a pseudo-grid to find pairs, based
* on the Bournaud technique
*/
#ifdef MULTIPHASE
int find_sticky_pairs(int ii)
{
int jj;
int ci,cj;
int i,j;
double dx,dy,dz,r2,r;
double dvx,dvy,dvz,vdotr,vdotr2,dv;
if (ii>=N_gas)
return -11;
ci = StickyIndex[ii].CellIndex;
if (ci>=Ncxyz)
return -12;
jj = ii;
do
{
jj++;
if (jj>=N_gas)
return -1;
cj = StickyIndex[jj].CellIndex;
if (cj>=Ncxyz)
return -2;
if (ci!=cj)
return -3;
if (SphP[StickyIndex[jj].Index].StickyFlag) /* if its free, return it, else, go to the next one */
{
/* here we check that we can interact */
i = StickyIndex[ii].Index;
j = StickyIndex[jj].Index;
dx = P[i].Pos[0] - P[j].Pos[0];
dy = P[i].Pos[1] - P[j].Pos[1];
dz = P[i].Pos[2] - P[j].Pos[2];
r2 = dx * dx + dy * dy + dz * dz;
r = sqrt(r2);
dvx = P[i].Vel[0] - P[j].Vel[0];
dvy = P[i].Vel[1] - P[j].Vel[1];
dvz = P[i].Vel[2] - P[j].Vel[2];
vdotr = dx * dvx + dy * dvy + dz * dvz;
vdotr2 = vdotr;
dv = vdotr2/r;
/*disable Mon Feb 21 16:47:28 CET 2011*/ //if ( (vdotr2<0) && (-vdotr2>All.StickyMinVelocity))
{
//printf("(1)%d %d\n",i,j);
return StickyIndex[jj].Index; /* ok, we accept the collision */
}
}
}
while(1);
}
#endif
void sticky_compute_energy_kin(int mode)
{
int i;
double DeltaEgyKin;
double Tot_DeltaEgyKin;
DeltaEgyKin = 0;
Tot_DeltaEgyKin = 0;
if (mode==1)
{
LocalSysState.EnergyKin1 = 0;
LocalSysState.EnergyKin2 = 0;
}
for(i = 0; i < N_gas; i++)
{
if (P[i].Type==0)
{
if (mode==1)
LocalSysState.EnergyKin1 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
else
LocalSysState.EnergyKin2 += 0.5 * P[i].Mass * (P[i].Vel[0]*P[i].Vel[0]+P[i].Vel[1]*P[i].Vel[1]+P[i].Vel[2]*P[i].Vel[2]);
}
}
if (mode==2)
{
DeltaEgyKin = LocalSysState.EnergyKin2 - LocalSysState.EnergyKin1;
MPI_Reduce(&DeltaEgyKin, &Tot_DeltaEgyKin, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
LocalSysState.EnergyRadSticky -= DeltaEgyKin;
}
}
/*! compute sticky collisions
*
*/
void sticky_collisions()
{
int i,j,k;
int startnode;
FLOAT *pos, *vel;
FLOAT mass, h_i;
int phase;
double vcm[3],newv1[3],newv2[3];
double dv1[3],dv2[3];
double dx, dy, dz, dvx, dvy, dvz;
double r,r2;
double vdotr, vdotr2;
double dv;
double beta_r,beta_t;
double M1M,M2M,M12;
double dbeta_dv,dbeta_dv_er_x,dbeta_dv_er_y,dbeta_dv_er_z;
double dv_beta_t_x,dv_beta_t_y,dv_beta_t_z;
NumColLocal = 0;
NumNoColLocal = 0;
beta_r = All.StickyBetaR;
beta_t = All.StickyBetaT;
double xi,yi,zi;
int ix,iy,iz;
size_t bytes;
if (ThisTask==0)
printf("Start sticky collisions...\n");
/*****************************/
/* first, update StickyIndex */
/*****************************/
if (All.StickyUseGridForCollisions)
{
/* allocate memory */
if(!(StickyIndex = malloc(bytes = N_gas * sizeof(struct Sticky_index))))
{
printf("failed to allocate memory for `StickyIndex' (%g MB).\n", bytes / (1024.0 * 1024.0));
endrun(2233);
}
/* loop gas particles */
for(i = 0; i < N_gas; i++)
{
//if (P[i].Type==0) /* we also take also into account the stars not transfered in stars now */
{
StickyIndex[i].Index = i;
StickyIndex[i].CellIndex = Ncxyz;
StickyIndex[i].Flag = 1;
if (SphP[i].StickyFlag) /* particle that may collide (set in phase.c) */
{
/* scale between 0 and nx,ny,nz */
xi = All.StickyGridNx*(P[i].Pos[0]-All.StickyGridXmin)/(All.StickyGridXmax-All.StickyGridXmin);
yi = All.StickyGridNy*(P[i].Pos[1]-All.StickyGridYmin)/(All.StickyGridYmax-All.StickyGridYmin);
zi = All.StickyGridNz*(P[i].Pos[2]-All.StickyGridZmin)/(All.StickyGridZmax-All.StickyGridZmin);
ix = (int)(xi);
iy = (int)(yi);
iz = (int)(zi);
if (ix >= 0 && ix < All.StickyGridNx)
if (iy >= 0 && iy < All.StickyGridNy)
if (iz >= 0 && iz < All.StickyGridNz)
{
j = (ix)*(Ncxy) + (iy)*(Ncz) + (iz);
StickyIndex[i].Index = i;
StickyIndex[i].CellIndex = j;
}
}
}
}
/* sort particles */
qsort(StickyIndex, N_gas, sizeof(struct Sticky_index), compare_sticky_index);
/* record index in SphP (not necessary, but simplifies) */
/* loop over cells */
for(i = 0; i < N_gas; i++)
{
//printf("%d %d\n",StickyIndex[i].CellIndex,StickyIndex[i].Index);
j = StickyIndex[i].Index;
SphP[j].StickyIndex = i;
}
}
/*****************************/
/* loop over all particles */
/*****************************/
for(i = 0; i < N_gas; i++)
{
//#ifdef SFR
// if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
//#else
// if(P[i].Ti_endstep == All.Ti_Current)
//#endif
#ifdef SFR
if(P[i].Type == 0)
#endif
{
if(SphP[i].Phase == GAS_STICKY)
{
if(SphP[i].StickyFlag)
{
pos = P[i].Pos;
vel = SphP[i].VelPred;
mass = P[i].Mass;
h_i = SphP[i].Hsml;
phase = SphP[i].Phase;
startnode = All.MaxPart;
if (All.StickyUseGridForCollisions)
j = find_sticky_pairs(SphP[i].StickyIndex);
else
j = ngb_treefind_sticky_collisions(&pos[0], h_i, phase, &startnode);
/* check that particles are in the same cell */
/*
if (j>=0 && j!=i)
{
int i1,i2;
xi = All.StickyGridNx*(P[i].Pos[0]-All.StickyGridXmin)/(All.StickyGridXmax-All.StickyGridXmin);
yi = All.StickyGridNy*(P[i].Pos[1]-All.StickyGridYmin)/(All.StickyGridYmax-All.StickyGridYmin);
zi = All.StickyGridNz*(P[i].Pos[2]-All.StickyGridZmin)/(All.StickyGridZmax-All.StickyGridZmin);
ix = (int)(xi);
iy = (int)(yi);
iz = (int)(zi);
i1 = (ix)*(Ncxy) + (iy)*(Ncz) + (iz);
xi = All.StickyGridNx*(P[i].Pos[0]-All.StickyGridXmin)/(All.StickyGridXmax-All.StickyGridXmin);
yi = All.StickyGridNy*(P[i].Pos[1]-All.StickyGridYmin)/(All.StickyGridYmax-All.StickyGridYmin);
zi = All.StickyGridNz*(P[i].Pos[2]-All.StickyGridZmin)/(All.StickyGridZmax-All.StickyGridZmin);
ix = (int)(xi);
iy = (int)(yi);
iz = (int)(zi);
i2 = (ix)*(Ncxy) + (iy)*(Ncz) + (iz);
if (i1!=i2)
{
printf("\n\nSomething odd here!\n");
endrun("999888777");
}
}
*/
if (j>=0 && j!=i) /* particles may collide */
{
if (j>=N_gas)
{
printf("sticky_collisions : j=%d,N_gas=%d j>N_gas !!!",j,N_gas);
exit(-1);
}
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;
r = sqrt(r2);
dvx = vel[0] - SphP[j].VelPred[0];
dvy = vel[1] - SphP[j].VelPred[1];
dvz = vel[2] - SphP[j].VelPred[2];
vdotr = dx * dvx + dy * dvy + dz * dvz;
vdotr2 = vdotr;
dv = vdotr2/r;
if (vdotr2<0) /* particles approaches !!! */
{
/* compute new velocities */
M12 = mass+P[j].Mass;
M1M = mass/M12;
M2M = P[j].Mass/M12;
/* compute velocity of the mass center */
for(k = 0; k < 3; k++)
vcm[k] = (mass*vel[k]+P[j].Mass*SphP[j].VelPred[k])/M12;
dbeta_dv = (beta_r-beta_t)*dv;
dbeta_dv_er_x = -dbeta_dv *dx/r;
dbeta_dv_er_y = -dbeta_dv *dy/r;
dbeta_dv_er_z = -dbeta_dv *dz/r;
dv_beta_t_x = dvx*beta_t;
dv_beta_t_y = dvy*beta_t;
dv_beta_t_z = dvz*beta_t;
newv1[0] = M2M * ( +dv_beta_t_x - dbeta_dv_er_x ) + vcm[0];
newv1[1] = M2M * ( +dv_beta_t_y - dbeta_dv_er_y ) + vcm[1];
newv1[2] = M2M * ( +dv_beta_t_z - dbeta_dv_er_z ) + vcm[2];
newv2[0] = M1M * ( -dv_beta_t_x + dbeta_dv_er_x ) + vcm[0];
newv2[1] = M1M * ( -dv_beta_t_y + dbeta_dv_er_y ) + vcm[1];
newv2[2] = M1M * ( -dv_beta_t_z + dbeta_dv_er_z ) + vcm[2];
/* new velocities */
for(k = 0; k < 3; k++)
{
dv1[k] = newv1[k] - SphP[i].VelPred[k];
dv2[k] = newv2[k] - SphP[j].VelPred[k];
SphP[i].VelPred[k] = newv1[k];
SphP[j].VelPred[k] = newv2[k];
P[i].Vel[k] += dv1[k];
P[j].Vel[k] += dv2[k];
}
/* set particles as non sticky-active */
SphP[i].StickyFlag = 0;
SphP[j].StickyFlag = 0;
SphP[i].StickyTime = All.Time + All.StickyIdleTime;
SphP[j].StickyTime = All.Time + All.StickyIdleTime;
/* record collision */
NumColLocal+=2;
}
}
}
}
}
}
/* write some statistics */
MPI_Allreduce(&NumColLocal, &NumCol, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
//MPI_Allreduce(&NumColPotLocal,&NumColPot, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&NumNoColLocal, &NumNoCol, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if(ThisTask == 0)
{
fprintf(FdSticky, "Step %d, Time: %g NumColPot: %d NumCol: %d NumNoCol: %d\n", All.NumCurrentTiStep, All.Time,(int)NumColPot,(int)NumCol,(int)NumNoCol);
fflush(FdSticky);
}
if (All.StickyUseGridForCollisions)
free(StickyIndex);
if (ThisTask==0)
printf("sticky collisions done.\n");
}
/*! This function is the driver routine for the calculation of hydrodynamical
* force and rate of change of entropy due to shock heating for all active
* particles .
*
* During the first loop, each particle compute j (a potentially colliding particle)
* During the second loop, it check if it can collide and collide if it can.
*
*/
void sticky_collisions2(int loop)
{
long long ntot, ntotleft;
int i, j, k, n, ngrp, maxfill, source, ndone;
int *nbuffer, *noffset, *nsend_local, *nsend, *numlist, *ndonelist;
int level, sendTask, recvTask, nexport, place;
double tstart, tend, sumt, sumcomm;
double timecomp = 0, timecommsumm = 0, timeimbalance = 0, sumimbalance;
MPI_Status status;
double dv[3];
#ifdef PERIODIC
boxSize = All.BoxSize;
boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
boxHalf_X = boxHalf * LONG_X;
boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
boxHalf_Y = boxHalf * LONG_Y;
boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
boxHalf_Z = boxHalf * LONG_Z;
boxSize_Z = boxSize * LONG_Z;
#endif
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
double v1m,v2m;
#endif
if(All.ComovingIntegrationOn)
{
/* Factors for comoving integration of hydro */
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);
hubble_a2 = All.Time * All.Time * hubble_a;
hubble_a_a = All.Time * hubble_a;
ainv = 1.0/All.Time;
fac_mu = pow(All.Time, 3 * (GAMMA - 1) / 2) / All.Time;
fac_egy = pow(All.Time, 3 * (GAMMA - 1));
fac_vsic_fix = hubble_a * pow(All.Time, 3 * GAMMA_MINUS1);
a3inv = 1 / (All.Time * All.Time * All.Time);
atime = All.Time;
#ifdef FEEDBACK
fac_pow = fac_egy*atime*atime;
#endif
}
else
{
hubble_a = hubble_a2 = hubble_a_a = ainv = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0;
#ifdef FEEDBACK
fac_pow = 1.0;
#endif
}
/* `NumSphUpdate' gives the number of particles on this processor that want a force update */
for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
{
#ifdef SFR
//if((P[n].Ti_endstep == All.Ti_Current) && (P[n].Type == 0))
#else
//if(P[n].Ti_endstep == All.Ti_Current)
#endif
if((SphP[n].Phase == GAS_STICKY) && (SphP[n].StickyFlag) )
{
if (loop==1)
{
SphP[n].StickyMaxFs = 0;
SphP[n].StickyMaxID =-1;
}
SphP[n].StickyNgb =-1;
NumSphUpdate++;
}
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
if (loop==1)
for(i = 0, All.TotN_stickyactive = 0; i < NTask; i++)
All.TotN_stickyactive += numlist[i];
free(numlist);
noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */
nbuffer = malloc(sizeof(int) * NTask);
nsend_local = malloc(sizeof(int) * NTask);
nsend = malloc(sizeof(int) * NTask * NTask);
ndonelist = malloc(sizeof(int) * NTask);
i = 0; /* first particle for this task */
ntotleft = ntot; /* particles left for all tasks together */
NumColLocal = 0;
NumNoColLocal = 0;
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeSticky - NTask; i++)
{
#ifdef SFR
//if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
//if(P[i].Ti_endstep == All.Ti_Current)
#endif
if((SphP[i].Phase == GAS_STICKY)&&(SphP[i].StickyFlag))
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
sticky_evaluate(i, 0,loop);
/* the particle is not exported if an ngb has been found locally */
if (SphP[i].StickyNgb!=-1)
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
for(k = 0; k < 3; k++)
{
StickyDataIn[nexport].Pos[k] = P[i].Pos[k];
StickyDataIn[nexport].Vel[k] = SphP[i].VelPred[k];
}
StickyDataIn[nexport].Mass = P[i].Mass;
StickyDataIn[nexport].Hsml = SphP[i].Hsml;
StickyDataIn[nexport].StickyMaxID = SphP[i].StickyMaxID;
StickyDataIn[nexport].StickyMaxFs = SphP[i].StickyMaxFs;
StickyDataIn[nexport].StickyNgb = SphP[i].StickyNgb;
StickyDataIn[nexport].ID = P[i].ID;
StickyDataIn[nexport].Index = i;
StickyDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
}
}
}
tend = second();
timecomp += timediff(tstart, tend);
qsort(StickyDataIn, nexport, sizeof(struct stickydata_in), sticky_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
tstart = second();
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeSticky)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* get the particles */
MPI_Sendrecv(&StickyDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct stickydata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&StickyDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct stickydata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
/* now do the imported particles */
tstart = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
sticky_evaluate(j, 1,loop);
tend = second();
timecomp += timediff(tstart, tend);
/* do a block to measure imbalance */
tstart = second();
MPI_Barrier(MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* get the result */
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeSticky)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
/* send the results */
MPI_Sendrecv(&StickyDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct stickydata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B,
&StickyDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct stickydata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B, MPI_COMM_WORLD, &status);
/* add the result to the particles */
for(j = 0; j < nsend_local[recvTask]; j++)
{
source = j + noffset[recvTask];
place = StickyDataIn[source].Index;
if (StickyDataPartialResult[source].StickyMaxFs>SphP[place].StickyMaxFs)
{
SphP[place].StickyMaxID = StickyDataPartialResult[source].StickyMaxID;
SphP[place].StickyMaxFs = StickyDataPartialResult[source].StickyMaxFs;
}
if (StickyDataPartialResult[source].StickyNgb!=-1) /* copy only if met a neighbor */
{
SphP[place].StickyNgb = StickyDataPartialResult[source].StickyNgb;
for(k = 0; k < 3; k++)
SphP[place].StickyNewVel[k]= StickyDataPartialResult[source].StickyNewVel[k];
}
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
level = ngrp - 1;
}
MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
for(j = 0; j < NTask; j++)
ntotleft -= ndonelist[j];
}
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
if (loop==2)
{
/* do final operations on results */
double norm_dv;
double max_dv;
int amax_dv;
max_dv = 0;
for(i = 0; i < N_gas; i++)
#ifdef SFR
//if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
//if(P[i].Ti_endstep == All.Ti_Current)
#endif
if((SphP[i].Phase == GAS_STICKY)&&(SphP[i].StickyFlag))
{
if (SphP[i].StickyNgb!=-1)
{
//printf("(%d) COLLISION : %d -> %d\n",ThisTask,P[i].ID,SphP[i].StickyNgb);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[0],SphP[i].StickyNewVel[0]);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[1],SphP[i].StickyNewVel[1]);
//printf(" : OldV %g NewV %g\n",SphP[i].VelPred[2],SphP[i].StickyNewVel[2]);
/* new velocities */
for(k = 0; k < 3; k++)
{
dv[k] = SphP[i].StickyNewVel[k] - SphP[i].VelPred[k];
SphP[i].VelPred[k] = SphP[i].StickyNewVel[k];
P[i].Vel[k] += dv[k];
}
norm_dv = sqrt( dv[0]*dv[0] + dv[1]*dv[1] + dv[2]*dv[2]);
if (norm_dv>max_dv)
{
max_dv = norm_dv;
amax_dv= P[i].ID;
}
/* set particles as non sticky-active */
SphP[i].StickyFlag = 0;
SphP[i].StickyTime = All.Time;
/* record collision */
NumColLocal++;
}
else
NumNoColLocal++;
}
printf("(%d) max_dv=%g amax_dv=%d\n",ThisTask,max_dv,amax_dv);
}
tend = second();
timecomp += timediff(tstart, tend);
/* write some statistics */
if(loop==2)
{
MPI_Allreduce(&NumColLocal, &NumCol, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&NumNoColLocal, &NumNoCol, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if(ThisTask == 0)
{
//fprintf(FdSticky, "Step %d, Time: %g NumColPot: %d NumCol: %d NumNoCol: %d\n", All.NumCurrentTiStep, All.Time,(int)NumColPot,(int)NumCol,(int)NumNoCol);
fprintf(FdSticky, "%10d %15g %10d %10d %10d %10d\n", All.NumCurrentTiStep, All.Time,(int)All.TotN_sticky,(int)All.TotN_stickyactive,(int)NumCol,(int)NumNoCol);
fflush(FdSticky);
}
}
if (ThisTask==0)
printf("sticky collisions done.\n");
}
/*!
*
* first sticky
*
*/
void sticky_evaluate(int target, int mode, int loop)
{
int j,k, n, startnode, numngb;
FLOAT *pos, *vel;
FLOAT mass;
FLOAT h_i;
double dx, dy, dz;
double dvx,dvy,dvz;
double vdotr, vdotr2,dv,dv2;
double vcm[3],newvel[3];
double h_i2;
double h_j, r, r2;
int phase=0;
int StickyMaxID;
int id;
double M1M,M2M,M12;
double dbeta_dv,dbeta_dv_er_x,dbeta_dv_er_y,dbeta_dv_er_z;
double dv_beta_t_x,dv_beta_t_y,dv_beta_t_z;
double beta_r,beta_t;
double iM12;
double W0,W1,W2,vdotr2r2;
double dvxadxa,dvyadya,dvzadza;
double vdotr2r2dx,vdotr2r2dy,vdotr2r2dz;
float fs,maxfs;
int id_max;
int id_ngb;
beta_r = All.StickyBetaR;
beta_t = All.StickyBetaT;
if(mode == 0)
{
pos = P[target].Pos;
vel = SphP[target].VelPred;
mass= P[target].Mass;
h_i = SphP[target].Hsml;
id = P[target].ID;
StickyMaxID = SphP[target].StickyMaxID;
maxfs = SphP[target].StickyMaxFs;
id_max= SphP[target].StickyMaxID;
id_ngb= SphP[target].StickyNgb;
}
else
{
pos = StickyDataGet[target].Pos;
vel = StickyDataGet[target].Vel;
mass= StickyDataGet[target].Mass;
h_i = StickyDataGet[target].Hsml;
id = StickyDataGet[target].ID;
StickyMaxID = StickyDataGet[target].StickyMaxID;
maxfs = StickyDataGet[target].StickyMaxFs;
id_max= StickyDataGet[target].StickyMaxID;
id_ngb=StickyDataGet[target].StickyNgb;
}
h_i2 = h_i * h_i;
for(k = 0; k < 3; k++)
newvel[k] = vel[k];
/* Now start sticky 1 for this particle */
startnode = All.MaxPart;
do
{
numngb = ngb_treefind_pairs(&pos[0], h_i, phase, &startnode);
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];
#ifdef PERIODIC /* find the closest image in the given box size */
if(dx > boxHalf_X)
dx -= boxSize_X;
if(dx < -boxHalf_X)
dx += boxSize_X;
if(dy > boxHalf_Y)
dy -= boxSize_Y;
if(dy < -boxHalf_Y)
dy += boxSize_Y;
if(dz > boxHalf_Z)
dz -= boxSize_Z;
if(dz < -boxHalf_Z)
dz += boxSize_Z;
#endif
r2 = dx * dx + dy * dy + dz * dz;
h_j = SphP[j].Hsml;
if(r2 < h_i2 || r2 < h_j * h_j)
{
r = sqrt(r2);
if(r > 0)
{
dvx = vel[0] - SphP[j].VelPred[0];
dvy = vel[1] - SphP[j].VelPred[1];
dvz = vel[2] - SphP[j].VelPred[2];
vdotr = dx * dvx + dy * dvy + dz * dvz;
if(All.ComovingIntegrationOn)
vdotr2 = vdotr + hubble_a2 * r2;
else
vdotr2 = vdotr;
dv = vdotr2/r;
if (loop==1)
{
if (vdotr2<0) /* particles approaches !!! */
{
/*
* the next lines are no longer needed if we do not select on
* the final relative velocities of particles
* /
// /* compute new velocities */
//
// M12 = mass+P[j].Mass;
// M1M = mass/M12;
// M2M = P[j].Mass/M12;
//
// /* compute velocity of the mass center */
// for(k = 0; k < 3; k++)
// vcm[k] = (mass*vel[k]+P[j].Mass*SphP[j].VelPred[k])/M12;
//
//
// dbeta_dv = (beta_r-beta_t)*dv;
// dbeta_dv_er_x = -dbeta_dv *dx/r;
// dbeta_dv_er_y = -dbeta_dv *dy/r;
// dbeta_dv_er_z = -dbeta_dv *dz/r;
//
// dv_beta_t_x = dvx*beta_t;
// dv_beta_t_y = dvy*beta_t;
// dv_beta_t_z = dvz*beta_t;
//
// newvel[0] = M2M * ( +dv_beta_t_x - dbeta_dv_er_x ) + vcm[0];
// newvel[1] = M2M * ( +dv_beta_t_y - dbeta_dv_er_y ) + vcm[1];
// newvel[2] = M2M * ( +dv_beta_t_z - dbeta_dv_er_z ) + vcm[2];
//
//
// /* check that dv is not too large */
// dv2 = (newvel[0]-vel[0])*(newvel[0]-vel[0])
// + (newvel[1]-vel[1])*(newvel[1]-vel[1])
// + (newvel[2]-vel[2])*(newvel[2]-vel[2]);
/* disable Mon Feb 21 16:55:22 CET 2011 */ //if (dv2<All.StickyMaxVelocity*All.StickyMaxVelocity)
{
fs = r/( (0.5*(h_i+h_j)+r) * (0.5*(h_i+h_j)+r));
if (fs>maxfs)
{
id_max = P[j].ID;
maxfs = fs;
}
}
}
}
else
{
if ((SphP[j].StickyMaxID == id) && (StickyMaxID == P[j].ID))
{
id_ngb = P[j].ID;
/* compute new velocities : new version : Pfenniger 2001.02.16*/
iM12 = 1.0/(mass+P[j].Mass);
vdotr2r2 = vdotr2 / r2;
dvxadxa = (-dvx*ainv) + (-dx*hubble_a_a);
dvyadya = (-dvy*ainv) + (-dy*hubble_a_a);
dvzadza = (-dvz*ainv) + (-dz*hubble_a_a);
W0 = iM12 * ( (1.0-beta_t)*(dvxadxa) + (beta_t-beta_r)*(vdotr2r2*(-dx)*ainv));
W1 = iM12 * ( (1.0-beta_t)*(dvyadya) + (beta_t-beta_r)*(vdotr2r2*(-dy)*ainv));
W2 = iM12 * ( (1.0-beta_t)*(dvzadza) + (beta_t-beta_r)*(vdotr2r2*(-dz)*ainv));
newvel[0] = vel[0] + P[j].Mass*W0*atime;
newvel[1] = vel[1] + P[j].Mass*W1*atime;
newvel[2] = vel[2] + P[j].Mass*W2*atime;
// /* compute new velocities : old standard version */
//
// M12 = mass+P[j].Mass;
// M1M = mass/M12;
// M2M = P[j].Mass/M12;
//
// /* compute velocity of the mass center */
// for(k = 0; k < 3; k++)
// vcm[k] = (mass*vel[k]+P[j].Mass*SphP[j].VelPred[k])/M12;
//
//
// dbeta_dv = (beta_r-beta_t)*dv;
// dbeta_dv_er_x = -dbeta_dv *dx/r;
// dbeta_dv_er_y = -dbeta_dv *dy/r;
// dbeta_dv_er_z = -dbeta_dv *dz/r;
//
// dv_beta_t_x = dvx*beta_t;
// dv_beta_t_y = dvy*beta_t;
// dv_beta_t_z = dvz*beta_t;
//
//
// newvel[0] = M2M * ( +dv_beta_t_x - dbeta_dv_er_x ) + vcm[0];
// newvel[1] = M2M * ( +dv_beta_t_y - dbeta_dv_er_y ) + vcm[1];
// newvel[2] = M2M * ( +dv_beta_t_z - dbeta_dv_er_z ) + vcm[2];
}
}
}
}
}
}
while(startnode >= 0);
/* Now collect the result at the right place */
if(mode == 0)
{
if(loop == 1)
{
SphP[target].StickyMaxID = id_max;
SphP[target].StickyMaxFs = maxfs;
}
SphP[target].StickyNgb = id_ngb;
for(k = 0; k < 3; k++)
SphP[target].StickyNewVel[k] = newvel[k];
}
else
{
if(loop == 1)
{
StickyDataResult[target].StickyMaxID = id_max;
StickyDataResult[target].StickyMaxFs= maxfs;
}
StickyDataResult[target].StickyNgb= id_ngb;
for(k = 0; k < 3; k++)
StickyDataResult[target].StickyNewVel[k] = newvel[k];
}
}
/*! This is a comparison kernel for a sort routine, which is used to group
* particles that are going to be exported to the same CPU.
*/
int sticky_compare_key(const void *a, const void *b)
{
if(((struct stickydata_in *) a)->Task < (((struct stickydata_in *) b)->Task))
return -1;
if(((struct stickydata_in *) a)->Task > (((struct stickydata_in *) b)->Task))
return +1;
return 0;
}
#endif
c4science · Help