Homec4science
Diffusion Gear (master) Browse

Gear/src/fof.c
1c0fbaf83a18master

Last Change
Raw File

fof.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <gsl/gsl_math.h>
#include "allvars.h"
#include "proto.h"
#ifdef FOF
#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
#define MAX_STARS_PER_PART 100
#ifdef PERIODIC
/*! Macro that maps a distance to the nearest periodic neighbour */
#define NEAREST(x) (((x)>boxhalf)?((x)-boxsize):(((x)<-boxhalf)?((x)+boxsize):(x)))
/*! Size of 3D lock-up table for Ewald correction force */
#define EN 64
/*! 3D lock-up table for Ewald correction to force and potential. Only one
* octant is stored, the rest constructed by using the symmetry
*/
static FLOAT fcorrx[EN + 1][EN + 1][EN + 1];
static FLOAT fcorry[EN + 1][EN + 1][EN + 1];
static FLOAT fcorrz[EN + 1][EN + 1][EN + 1];
static FLOAT potcorr[EN + 1][EN + 1][EN + 1];
static double fac_intp;
#endif
static int NpHead;
static int NHead;
static double a3,a3inv,hubble_a,DeltaTime;
static struct fof_phdata
{
int index;
int prev;
int cpuprev;
FLOAT density;
}
*fofd;
static struct fof_topgroups_exchange
{
int Head; /*!< head index */
int Task; /*!< task hosting the true head */
}
*TopTSfrGroups,*TopTSfrGroups_local;
void fof(void)
{
double dt,dt_Myr;
int do_FOF;
int Ti1,Ti2;
/* check if it is time to do an fof */
if (All.ComovingIntegrationOn)
{
Ti1 = log(All.FoF_TimeLastFoF/All.TimeBegin) / All.Timebase_interval;
Ti2 = All.Ti_Current;
dt = get_cosmictime_difference(Ti1,Ti2);
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);
a3 = All.Time * All.Time * All.Time;
a3inv = 1/a3;
}
else
{
dt = All.Time-All.FoF_TimeLastFoF;
hubble_a = a3 = a3inv = 1;
}
if (ThisTask==0)
{
dt_Myr = dt/3.1536e13*All.UnitTime_in_s/All.HubbleParam;
printf("FoF: Time=%g LastTime=%g dt=%g (%g Myr, hubble_a=%g)\n",All.Time,All.FoF_TimeLastFoF,dt,dt_Myr,hubble_a);
}
DeltaTime = dt;
//if((All.Time - All.FoF_TimeLastFoF) >= All.FoF_TimeBetFoF)
if(dt>=All.FoF_TimeBetFoF)
{
if(ThisTask == 0)
{
printf("FoF: Starting...\n");
printf("FoF: last FoF=%g\n",All.FoF_TimeLastFoF);
fflush(stdout);
}
fof_init();
//fof_write_particles_id_and_indicies();
fof_find_densest_neighbour();
fof_find_local_heads();
fof_link_particles_to_local_head();
fof_clean_local_groups();
fof_compute_local_tails();
//printf("(%d) NHead=%08d NpHead=%08d\n",ThisTask,NHead,NpHead);
fof_regroup_pseudo_heads();
//fof_group_info();
//printf("(%d) NHead=%08d NpHead=%08d\n",ThisTask,NHead,NpHead);
//fof_write_local_groups();
fof_send_and_link_pseudo_heads();
fof_regroup_exported_pseudo_heads();
//fof_write_all_groups();
fof_allocate_groups();
fof_compute_groups_properties();
#ifdef SFR
//fof_star_formation();
fof_star_formation_and_IMF_sampling();
#endif
//fof_write_groups_properties();
fof_free_groups();
All.FoF_TimeLastFoF = All.Time;
if(ThisTask == 0)
{
printf("FoF: next FoF=%g \n",All.FoF_TimeLastFoF+All.FoF_TimeBetFoF);
printf("FoF: done.\n");
fflush(stdout);
}
}
}
void fof_init(void)
{
#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
}
void fof_find_densest_neighbour(void)
{
if(ThisTask == 0)
{
printf("FoF: find densest neighbour\n");
fflush(stdout);
}
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;
MPI_Status status;
/*
for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
{
if(P[n].Type == 0)
{
SphP[n].FOF_Head=n;
SphP[n].FOF_DensMax=-1;
SphP[n].FOF_CPUHead=-1;
}
}
N_gas=10;
*/
//compute_potential();
/* `NumSphUpdate' gives the number of particles on this processor that want a force update */
for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
{
if ( (P[n].Type == 0) && (SphP[n].Density*a3inv>All.FoF_ThresholdDensity) )
{
NumSphUpdate++;
SphP[n].FOF_Len=1;
SphP[n].FOF_Head=-1;
SphP[n].FOF_Tail=-1;
SphP[n].FOF_Next=-1;
SphP[n].FOF_Prev=-1;
SphP[n].FOF_CPUHead=-1;
SphP[n].FOF_CPUTail=-1;
SphP[n].FOF_CPUNext=-1;
SphP[n].FOF_CPUPrev=-1;
SphP[n].FOF_DensMax=-1;
SphP[n].FOF_Done=0;
}
}
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];
free(numlist);
if (ThisTask==0)
printf("FoF: Number of gas particles above threshold : %d\n",ntot);
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 */
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeFOF - NTask; i++)
if((P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity))
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
fof_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
for(k = 0; k < 3; k++)
{
FOFDataIn[nexport].Pos[k] = P[i].Pos[k];
}
FOFDataIn[nexport].Hsml = SphP[i].Hsml;
FOFDataIn[nexport].FOF_Prev = SphP[i].FOF_Prev;
FOFDataIn[nexport].FOF_DensMax = SphP[i].FOF_DensMax;
FOFDataIn[nexport].FOF_CPUPrev = SphP[i].FOF_CPUPrev;
FOFDataIn[nexport].Index = i;
FOFDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
}
}
qsort(FOFDataIn, nexport, sizeof(struct FOFdata_in), fof_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
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.BunchSizeFOF)
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(&FOFDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFdata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&FOFDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFdata_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];
}
/* now do the imported particles */
for(j = 0; j < nbuffer[ThisTask]; j++)
fof_evaluate(j, 1);
/* do a block to measure imbalance */
MPI_Barrier(MPI_COMM_WORLD);
/* get the result */
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.BunchSizeFOF)
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(&FOFDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFdata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B,
&FOFDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFdata_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 = FOFDataIn[source].Index;
if (FOFDataPartialResult[source].FOF_DensMax>SphP[place].FOF_DensMax)
{
SphP[place].FOF_DensMax = FOFDataPartialResult[source].FOF_DensMax;
SphP[place].FOF_Prev = FOFDataPartialResult[source].FOF_Prev;
SphP[place].FOF_CPUPrev = FOFDataPartialResult[source].FOF_CPUPrev;
}
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
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);
}
/*! This function is the 'core' of the FOF computation.
*
* For each particle, we find its densest neighboring particle
* and store its index and cpu where its reside in the variables
*
* FOF_Prev = j;
* FOF_CPUPrev = ThisTask;
*
*/
void fof_evaluate(int target, int mode)
{
int j, n, startnode, numngb;
FLOAT *pos;
FLOAT h_i;
int phase=0;
double dx, dy, dz;
double r2, h2;
int FOF_Prev;
FLOAT FOF_DensMax;
int FOF_CPUPrev;
int id;
if(mode == 0)
{
pos = P[target].Pos;
h_i = SphP[target].Hsml;
FOF_Prev = SphP[target].FOF_Prev;
FOF_DensMax = SphP[target].FOF_DensMax;
FOF_CPUPrev = SphP[target].FOF_CPUPrev;
id = P[target].ID;
}
else
{
pos = FOFDataGet[target].Pos;
h_i = FOFDataGet[target].Hsml;
FOF_Prev = FOFDataGet[target].FOF_Prev;
FOF_DensMax = FOFDataGet[target].FOF_DensMax;
FOF_CPUPrev = FOFDataGet[target].FOF_CPUPrev;
id = FOFDataGet[target].ID;
}
h_i = h_i*All.FoF_HsmlSearchRadiusFactor;
h2 = h_i*h_i;
/* Now start the actual SPH computation for this particle */
startnode = All.MaxPart;
do
{
numngb = ngb_treefind_variable(&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;
if (r2 < h2)
{
if (SphP[j].Density>FOF_DensMax)
{
FOF_DensMax = SphP[j].Density;
FOF_Prev = j;
FOF_CPUPrev = ThisTask;
}
}
}
}
while(startnode >= 0);
/* Now collect the result at the right place */
if(mode == 0)
{
SphP[target].FOF_DensMax = FOF_DensMax;
SphP[target].FOF_Prev = FOF_Prev;
SphP[target].FOF_CPUPrev = FOF_CPUPrev;
}
else
{
FOFDataResult[target].FOF_DensMax = FOF_DensMax;
FOFDataResult[target].FOF_Prev = FOF_Prev;
FOFDataResult[target].FOF_CPUPrev = FOF_CPUPrev;
}
}
/*! Compute the number of true head and pseudo head
*
* A true head is a particle for which the densest neighboring particle
* is itself (residing in the same cpu).
*
* (SphP[i].FOF_Prev==i) && (SphP[i].FOF_CPUPrev==ThisTask)
*
* -> SphP[i].FOF_Head = i;
* -> SphP[i].FOF_CPUHead = ThisTask;
*
* A pseudo head is a particle for which the densest neighboring particle
* is on another CPU
*
* SphP[i].FOF_CPUPrev!=ThisTask
*
* -> SphP[i].FOF_Head = i;
* -> SphP[i].FOF_CPUHead = SphP[i].FOF_CPUPrev;
*
*
* New definitions FROM HERE:
*
* a head is : SphP[i].FOF_Head==i
* a true head is : (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead==ThisTask)
* a pseudo head is : (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead!=ThisTask)
*
*/
void fof_find_local_heads(void)
{
int i;
NHead=0;
NpHead=0;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
if ( (SphP[i].FOF_Prev==i) && (SphP[i].FOF_CPUPrev==ThisTask) ) /* a true head */
{
SphP[i].FOF_Head = i;
SphP[i].FOF_CPUHead = ThisTask;
NHead++;
}
if ( (SphP[i].FOF_CPUPrev!=ThisTask) && (SphP[i].FOF_CPUPrev>=0) ) /* if paticles needs to be exported */
{
SphP[i].FOF_Head = i; /* assume these particles to be (local) head */
SphP[i].FOF_CPUHead = SphP[i].FOF_CPUPrev;
NpHead++;
}
}
else
{
SphP[i].FOF_Head = -1;
SphP[i].FOF_Prev = -1;
}
}
/*! Link each particle to a head
*
*
*/
void fof_link_particles_to_local_head(void)
{
int i;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
fof_link_to_local_head(i);
}
}
/*! Link each particle to a head
*
* Each particle i try to find recursively
* a particle that is a head.
*
* The variable SphP[i].FOF_Next is used here
* to link particles with their descendents.
*
*/
void fof_link_to_local_head(int i)
{
/*
FOF_Prev==-1 : dead end
FOF_Prev==i : head particle
FOF_Head==i : head particle
FOF_Next==-1 : particle at the queue
*/
int k;
int n;
double dx,dy,dz,r2;
if (SphP[i].FOF_Head!=-1) /* the particle is already in a list */
return;
if (SphP[i].FOF_Prev==-1) /* mark it as a dead end */
return;
/* some physical restrictions */
if (SphP[i].Density*a3inv<All.FoF_ThresholdDensity)
{
SphP[i].FOF_Prev==-1; /* mark it as a dead end */
SphP[i].FOF_Head==-1;
return;
}
/* now, look for previous particle */
k = SphP[i].FOF_Prev;
/* the previous particle is ok, but is not included in a list */
if (SphP[k].FOF_Head==-1)
fof_link_to_local_head(k); /* find a head for this particle */
/* at this point k has a head or is a dead end */
if (SphP[k].FOF_Prev==-1) /* dead end */
{
SphP[i].FOF_Prev=-1;
SphP[i].FOF_Head=-1;
return;
}
else /* we can include it to the queue */
{
if (SphP[k].FOF_Head==-1)
endrun(1953);
/* !!! particle-particle comparison is not allowed here, as that may not be respected with multiple cpu !!! */
// if(All.FoF_HsmlMaxDistFactor>=0)
// {
//
// /* check the distance */
//
// dx = P[i].Pos[0] - P[k].Pos[0];
// dy = P[i].Pos[1] - P[k].Pos[1];
// dz = P[i].Pos[2] - P[k].Pos[2];
//
//
//#ifdef PERIODIC
// 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;
//
// if (r2 > pow(All.FoF_HsmlMaxDistFactor* SphP[i].Hsml,2)) /* the particle is too far away */
// {
// SphP[i].FOF_Prev=-1; /* mark it as a dead end */
// SphP[i].FOF_Head=-1;
// return;
// }
// }
/* now, everything is ok, add it to the list */
SphP[i].FOF_Head = SphP[k].FOF_Head;
/* deal with next */
n = SphP[k].FOF_Next;
SphP[k].FOF_Next = i;
SphP[i].FOF_Prev = k;
if (n!=-1) /* correct links with next particle */
{
SphP[i].FOF_Next = n;
SphP[n].FOF_Prev = i;
}
}
}
/*! Remove true heads based on some physical critteria
*
* Particles from groups that does not fit some physical critteria
* are marked as dead end, i.e.:
*
* SphP[next].FOF_Prev=-1;
* SphP[next].FOF_Head=-1;
*
*/
void fof_clean_local_groups(void)
{
int i;
int next;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if ( (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead==ThisTask) ) /* a true head */
{
//printf("a head : i=%d id=%08d (%g %g %g) rho=%g next=%d\n",i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],SphP[i].Density,SphP[i].FOF_Next);
/* here, we could select a group according to its central density or the number of particles linked to it */
if (SphP[i].Density*a3inv<All.FoF_MinHeadDensity)
{
SphP[i].FOF_Prev=-1; /* flag as dead end */
SphP[i].FOF_Head=-1;
NHead--;
next=i;
while(next!=-1)
{
SphP[next].FOF_Prev=-1; /* flag as dead end */
SphP[next].FOF_Head=-1;
next=SphP[next].FOF_Next;
}
}
}
}
/*! Compute local tails
*
* For each head (true or pseudo), loop over particles from the groups
* and record the tail into the variable :
*
* SphP[i].FOF_Tail
*
*/
void fof_compute_local_tails(void)
{
int i;
int tail;
int next;
FLOAT Density;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if ( (SphP[i].FOF_Head==i) ) /* a head or pseudo head */
{
next=i;
Density = SphP[next].Density;
while(next!=-1)
{
if (SphP[next].Density>Density)
{
printf("(%d) i=%d iID=%d next=%d nextID=%d Density=%g SphP[next].Density=%g\n",ThisTask,i,P[i].ID,next,P[next].ID,Density,SphP[next].Density);
endrun(1255);
}
tail = next;
next=SphP[next].FOF_Next;
}
SphP[i].FOF_Tail = tail;
}
}
/*! Regroups pseudo head
*
* Pseudo heads that points towards the same distant particle
* are grouped
*
* SphP[i].FOF_Tail
*
*/
void fof_regroup_pseudo_heads(void)
{
int i,it;
int head,prev,cpuprev,tail,next;
int OldNpHead;
if (ThisTask==0)
printf("FoF: regroup pseudo heads\n");
printf("FoF: (%d) NpHead = %d\n",ThisTask,NpHead);
if (NpHead>0)
{
fofd = malloc(sizeof(struct fof_phdata) * NpHead);
for(i = 0,it = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
if ((SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead!=ThisTask)) /* a pseudo head */
{
fofd[it].index = i;
fofd[it].prev = SphP[i].FOF_Prev;
fofd[it].cpuprev = SphP[i].FOF_CPUPrev;
fofd[it].density = SphP[i].Density;
it++;
}
}
qsort(fofd, NpHead, sizeof(struct fof_phdata), fof_compare_phdata);
/* now, loop over the sorted structure */
head = fofd[0].index;
prev = fofd[0].prev;
cpuprev = fofd[0].cpuprev;
OldNpHead = NpHead;
for(it = 1; it < OldNpHead; it++)
{
if( (fofd[it].prev != prev) || (fofd[it].cpuprev != cpuprev)) /* the pseudo head points towards another particle */
{
head = fofd[it].index; /* it becomes the new head */
prev = fofd[it].prev;
cpuprev = fofd[it].cpuprev;
}
else
{
NpHead--;
i = fofd[it].index;
if ( SphP[i].Density > SphP[head].Density) /* this allow the particle to be in density order */
{
/* the particle becomes the head */
i = head;
head = fofd[it].index;
prev = fofd[it].prev;
}
/* we can add the particle (i) to the tail of the pseudo head (head) */
/* change the head */
next = i;
while(next!=-1)
{
SphP[next].FOF_Head = head;
next = SphP[next].FOF_Next;
}
tail = SphP[head].FOF_Tail; /* !!! here, the density order is no longer conserved !!! */
SphP[head].FOF_Tail = SphP[i].FOF_Tail;
/* deal with tail */
SphP[tail].FOF_Next = i;
SphP[i].FOF_Prev = tail;
}
}
free(fofd);
}
}
/*! Regroups exported pseudo head
*
* 1) Link Pseudo heads (that have been exported) to local groups if needed
* 2) Pseudo heads (that have been exported) that points towards the same distant head
* are grouped
*
* SphP[i].FOF_Tail
*
*/
void fof_regroup_exported_pseudo_heads(void)
{
int i,it;
int head,prev,cpuprev,tail,next;
int OldNpHead;
if (ThisTask==0)
printf("FoF: regroup exported pseudo heads\n");
for(i = 0,it = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
if ((SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUPrev==-2)) /* a true head, but among them, some may be old spseudo head */
{
if ( (SphP[i].FOF_CPUHead==ThisTask) ) /* its new head is now local, so, link it to local head */
{
NpHead--;
/* link it to the local true head */
head = SphP[i].FOF_Prev; /* as it is still a pseudo head */
/* change the head */
next = i;
while(next!=-1)
{
SphP[next].FOF_Head = head;
next = SphP[next].FOF_Next;
}
tail = SphP[head].FOF_Tail; /* !!! here, the density order is no longer conserved !!! */
SphP[head].FOF_Tail = SphP[i].FOF_Tail;
/* deal with tail */
SphP[tail].FOF_Next = i;
SphP[i].FOF_Prev = tail;
}
else
{
SphP[i].FOF_CPUPrev = SphP[i].FOF_CPUHead;
}
}
}
if (NpHead>0)
{
fofd = malloc(sizeof(struct fof_phdata) * NpHead);
NpHead = 0; /* re-count as some pHead have change their status after exportation */
for(i = 0,it = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
if ((SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead!=ThisTask)) /* a pseudo head */
{
fofd[it].index = i;
fofd[it].prev = SphP[i].FOF_Prev;
fofd[it].cpuprev = SphP[i].FOF_CPUPrev;
fofd[it].density = SphP[i].Density;
it++;
NpHead++;
}
}
qsort(fofd, NpHead, sizeof(struct fof_phdata), fof_compare_phdata);
/* now, loop over the sorted structure */
head = fofd[0].index;
prev = fofd[0].prev;
cpuprev = fofd[0].cpuprev;
OldNpHead = NpHead;
for(it = 1; it < OldNpHead; it++)
{
if( (fofd[it].prev != prev) || (fofd[it].cpuprev != cpuprev)) /* the pseudo head points towards another particle */
{
head = fofd[it].index; /* it becomes the new head */
prev = fofd[it].prev;
cpuprev = fofd[it].cpuprev;
}
else /* the particle share the same head */
{
NpHead--;
i = fofd[it].index;
if ( SphP[i].Density > SphP[head].Density) /* this allow the particle to be in density order */
{
/* the particle becomes the head */
i = head;
head = fofd[it].index;
prev = fofd[it].prev;
}
/* we can add the particle (i) to the tail of the pseudo head (head) */
/* change the head */
next = i;
while(next!=-1)
{
SphP[next].FOF_Head = head;
next = SphP[next].FOF_Next;
}
tail = SphP[head].FOF_Tail; /* !!! here, the density order is no longer conserved !!! */
SphP[head].FOF_Tail = SphP[i].FOF_Tail;
/* deal with tail */
SphP[tail].FOF_Next = i;
SphP[i].FOF_Prev = tail;
}
}
free(fofd);
}
}
void fof_write_particles_id_and_indicies(void)
{
FILE *fd;
char fname[500];
int i;
sprintf(fname, "%d_snap.lst",ThisTask);
fd = fopen(fname, "w");
printf(".... writing particles...\n");
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
{
fprintf(fd,"%08d %08d\n",P[i].ID,i);
}
fclose(fd);
}
void fof_write_local_groups(void)
{
int i;
FILE *fd;
char fname[500];
int next;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if ( (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead==ThisTask) ) /* a true head */
{
printf("(%d) a head : i=%d id=%08d (%g %g %g) rho=%g next=%d\n",ThisTask,i,P[i].ID,P[i].Pos[0],P[i].Pos[1],P[i].Pos[2],SphP[i].Density,SphP[i].FOF_Next);
sprintf(fname, "groups/%d_group%08d.lst",ThisTask,P[i].ID);
fd = fopen(fname, "w");
next=i;
while(next!=-1)
{
fprintf(fd,"%08d\n",P[next].ID);
next=SphP[next].FOF_Next;
}
fclose(fd);
}
}
void fof_write_all_groups(void)
{
int i;
FILE *fd;
char fname[500];
int next;
if (ThisTask==0)
printf("FoF: write all groups\n");
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if (SphP[i].FOF_Head==i)
{
if (SphP[i].FOF_CPUHead==ThisTask) /* a true head */
{
//printf("(%d) a head : i=%d head=%d on %d \n",ThisTask,i,SphP[i].FOF_Head,SphP[i].FOF_CPUHead);
sprintf(fname, "groups/%08dat%02d_group%08dat%02d.lst",i,ThisTask,SphP[i].FOF_Head,SphP[i].FOF_CPUHead);
}
else
{ /* a pseudo head */
//printf("(%d) a head (p) : i=%d head=%d on %d \n",ThisTask,i,SphP[i].FOF_Prev,SphP[i].FOF_CPUHead,SphP[i].FOF_CPUHead);
sprintf(fname, "groups/%08dat%02d_group%08dat%02d.lst",i,ThisTask,SphP[i].FOF_Prev,SphP[i].FOF_CPUHead);
}
fd = fopen(fname, "w");
next=i;
while(next!=-1)
{
fprintf(fd,"%08d\n",P[next].ID);
next=SphP[next].FOF_Next;
}
fclose(fd);
}
}
void fof_compute_groups_properties(void)
{
int i,j,k;
int ng;
int next;
int level,ngrp;
int sendTask,recvTask;
int *nsend_local, *nsend,*noffset;
int place;
int maxig,gidx;
int ntg, npg;
MPI_Status status;
if (ThisTask==0)
printf("FoF: compute groups properties\n");
noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */
nsend_local = malloc(sizeof(int) * NTask);
nsend = malloc(sizeof(int) * NTask * NTask);
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
ng=0;
ntg = 0;
npg = 0;
/* loop over heads and create the groups */
for(i = 0; i < N_gas; i++)
{
SphP[i].FOF_gidx=-1; /* local indice of the group */
SphP[i].FOF_gid=-1; /* id of the group */
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if (SphP[i].FOF_Head==i) /* the particle is a head */
{
if (ng==(NpHead+NHead))
{
printf("(%d) FoF: warning g=%d >= %d+%d !\n",ThisTask,ng,NpHead,NHead);
endrun(1961);
}
if (SphP[i].FOF_CPUHead==ThisTask) /* a true head */
{
//printf("(%d) a head : head i=%d on %d \n",ThisTask,SphP[i].FOF_Head,SphP[i].FOF_CPUHead);
Groups[ng].Head = SphP[i].FOF_Head;
Groups[ng].HeadID = P[i].ID;
for (k=0;k<3;k++)
Groups[ng].HeadPos[k] = P[i].Pos[k];
ntg++;
}
else
{ /* a pseudo head */
//printf("(%d) a head (p) : head i=%d on %d \n",ThisTask,SphP[i].FOF_Prev,SphP[i].FOF_CPUHead,SphP[i].FOF_CPUHead);
Groups[ng].Head = SphP[i].FOF_Prev;
Groups[ng].HeadID = -1; /* don't know the ID */
npg++;
}
Groups[ng].Task = SphP[i].FOF_CPUHead;
Groups[ng].LocalHead = i;
nsend_local[SphP[i].FOF_CPUHead]++; /* count number of groups exported towards proc FOF_CPUHead */
/* some init */
fof_init_group_properties(ng);
ng++;
}
}
Ngroups = ng;
Ntgroups = ntg;
Npgroups = npg;
MPI_Allreduce(&Ngroups, &Tot_Ngroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&Ntgroups, &Tot_Ntgroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&Npgroups, &Tot_Npgroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (ThisTask==0)
{
printf("FoF: %d groups selected for star formation.\n",Tot_Ngroups);
printf("FoF: %d true groups. \n",Tot_Ntgroups);
printf("FoF: %d pseudo groups. \n",Tot_Npgroups);
}
/* sort groups */
qsort(Groups, Ngroups, sizeof(struct fof_groups_data), fof_compare_groups_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
/* gather number of exported groups */
/* info : nsend[j * NTask + i] : number of elements of j send to i */
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
/* find the max number of imported groups */
for (i=0,maxig=0;i<NTask;i++)
{
for (j=0;j<NTask;j++)
{
/* number of groups received by i */
if (i!=j)
maxig = dmax(maxig ,nsend[j * NTask + i] );
}
}
GroupsGet = malloc(sizeof(struct fof_groups_data) * maxig);
GroupsRecv = malloc(sizeof(struct fof_groups_data) * maxig);
/* link true head with their group */
for (i=0;i<ng;i++)
{
Groups[i].Index = i; /* compute index */
if( Groups[i].Task == ThisTask ) /* a true head */
SphP[Groups[i].Head].FOF_gidx = i; /* link towards the group id */ /*this is no longer needed, as done after */
}
/*
/* first communication (send peudo groups)
/*
/* get Head ID and Head Position
/*
*/
for(level = 1; level < (1 << PTask); level++)
{
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if((recvTask < NTask))
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
if (nsend[recvTask * NTask + ThisTask]>maxig)
endrun(1963);
/* get the groups */
MPI_Sendrecv(&Groups[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&GroupsGet[0],
nsend[recvTask * NTask + ThisTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status);
/* do the computation */
for(j = 0; j < nsend[recvTask * NTask + ThisTask]; j++)
{
if (GroupsGet[j].Task!=ThisTask)
endrun(1965);
GroupsGet[j].HeadID = P[GroupsGet[j].Head].ID;
for (k=0;k<3;k++)
GroupsGet[j].HeadPos[k] = P[GroupsGet[j].Head].Pos[k];
}
/* get the result */
MPI_Sendrecv(&GroupsGet[0],
nsend[recvTask * NTask + ThisTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&GroupsRecv[0],
nsend_local[recvTask] * sizeof(struct fof_groups_data),MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status
);
/* add the result to the groups */
for(j = 0; j < nsend_local[recvTask]; j++)
{
place = GroupsRecv[j].Index;
Groups[place].HeadID = GroupsRecv[j].HeadID;
for (k=0;k<3;k++)
Groups[place].HeadPos[k] = GroupsRecv[j].HeadPos[k];
}
}
}
}
level = ngrp - 1;
}
/*
/* loop over all particles in a group and attribute the group ID (this is needed for the groups potential computation)
/* The group ID correspond to the ID of the head particle.
*/
for (i=0;i<Ngroups;i++)
{
next=Groups[i].LocalHead; /* indice towards Head or Pseudo Head */
while(next!=-1) /* loop over all local members */
{
SphP[next].FOF_gidx=i; /* link all particles to its local group (true or pseudo) */
SphP[next].FOF_gid = Groups[i].HeadID; /* set the group ID */
next=SphP[next].FOF_Next;
}
}
#ifndef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
/* now here, we can compute potential, in order to estimate the virial fraction for each group */
if (Tot_Ngroups>0)
{
if (ThisTask==0)
printf("FoF: ");
compute_potential();
fof_compute_groups_potential();
}
#endif
/*
/* compute properties of groups using local particles
*/
for (i=0;i<Ngroups;i++)
{
next=Groups[i].LocalHead; /* indice towards Head or Pseudo Head */
Groups[i].Nlocal=0;
while(next!=-1) /* loop over all local members */
{
fof_add_particle_properties(next,i,0);
/* count local particles linked to the head */
Groups[i].Nlocal++;
next=SphP[next].FOF_Next;
}
}
/*
/* second communication (send peudo groups)
*/
/* now, need to communicate */
for(level = 1; level < (1 << PTask); level++)
{
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if((recvTask < NTask))
{
//printf("(%02d) -> %02d (level=%02d ngrp=%02d)\n",sendTask,recvTask,level,ngrp);
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
if (nsend[recvTask * NTask + ThisTask]>maxig)
endrun(1963);
/* get the particles */
MPI_Sendrecv(&Groups[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&GroupsGet[0],
nsend[recvTask * NTask + ThisTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status);
for(j = 0; j < nsend[recvTask * NTask + ThisTask]; j++)
{
//if (ThisTask==0)
// printf("(%d) head=%08d \n",ThisTask,GroupsGet[j].Head);
if (GroupsGet[j].Task!=ThisTask)
endrun(1967);
gidx = SphP[ GroupsGet[j].Head ].FOF_gidx ;
if (Groups[gidx].Head!=GroupsGet[j].Head)
{
printf("(%d) %d -> %d\n ",ThisTask,sendTask,recvTask);
printf("(%d) something is wrong here... GroupsGet[j].Head=%d GroupsGet[j].Task=%d gid=%d Groups[gid].Head=%d\n",ThisTask, GroupsGet[j].Head,GroupsGet[j].Task,gidx,Groups[gidx].Head);
endrun(1968);
}
fof_add_particle_properties(j,gidx,1);
}
}
}
}
level = ngrp - 1;
}
fof_final_operations_on_groups_properties();
fof_set_groups_for_sfr();
/* now, each pseudo groups needs to be exported */
/* WARNING : in term of communication, this part is not optimized at all
as we communicte the whole structure while only some items are needed
*/
for(level = 1; level < (1 << PTask); level++)
{
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if((recvTask < NTask))
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
if (nsend[recvTask * NTask + ThisTask]>maxig)
endrun(1963);
/* get the groups */
MPI_Sendrecv(&Groups[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&GroupsGet[0],
nsend[recvTask * NTask + ThisTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status);
/* do the computation */
for(j = 0; j < nsend[recvTask * NTask + ThisTask]; j++)
{
if (GroupsGet[j].Task!=ThisTask)
endrun(1965);
GroupsGet[j].SfrFlag = Groups[ SphP[ GroupsGet[j].Head ].FOF_gidx ].SfrFlag;
}
/* get the result */
MPI_Sendrecv(&GroupsGet[0],
nsend[recvTask * NTask + ThisTask] * sizeof(struct fof_groups_data), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&GroupsRecv[0],
nsend_local[recvTask] * sizeof(struct fof_groups_data),MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status
);
/* add the result to the groups */
for(j = 0; j < nsend_local[recvTask]; j++)
{
place = GroupsRecv[j].Index;
Groups[place].SfrFlag = GroupsRecv[j].SfrFlag;
}
}
}
}
level = ngrp - 1;
}
/* count star forming groups */
Nsfrgroups=0;
if (Ngroups>0)
{
for (gidx=0;gidx<Ngroups;gidx++)
{
if( Groups[gidx].SfrFlag == 1 ) /* the group must be flagged as forming stars */
Nsfrgroups++;
}
}
MPI_Allreduce(&Nsfrgroups, &Tot_Nsfrgroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (ThisTask==0)
printf("FoF: %d groups flagged for star formation.\n",Tot_Nsfrgroups);
free(nsend);
free(nsend_local);
free(noffset);
free(GroupsRecv);
free(GroupsGet);
}
void fof_init_group_properties(int gid)
{
int k;
Groups[gid].Mass = 0;
for (k=0;k<3;k++)
{
Groups[gid].MassCenter[k] = 0;
Groups[gid].MV2[k] = 0;
Groups[gid].MV[k] = 0;
}
Groups[gid].DensityMax = 0;
Groups[gid].DensityMin = pow(2,32);
Groups[gid].PotentialMax = 0;
Groups[gid].RadiusMax = 0;
Groups[gid].W = 0;
Groups[gid].K = 0;
Groups[gid].U = 0;
Groups[gid].N = 0;
}
/*! Include particle properties into a group
* if mode=0, we add the properties of a local particle
* if mode=1, we add the properties of a group (containing the contribution of the local particles)
*
*/
void fof_add_particle_properties(int i, int gid, int mode)
{
int k;
int N;
float Mass;
float DensityMax;
float DensityMin;
float PotentialMax;
float MassCenter[3];
float MV[3];
float MV2[3];
float W;
float U;
float RadiusMax;
float dx,dy,dz;
if (mode==0)
{
N = 1;
Mass = P[i].Mass;
DensityMax = SphP[i].Density;
DensityMin = SphP[i].Density;
PotentialMax = P[i].Mass*P[i].Potential;
W = P[i].Mass*P[i].Potential;
#ifdef ISOTHERM_EQS
U = P[i].Mass*SphP[i].Entropy;
#else
#ifdef MULTIPHASE
if (SphP[i].Phase== GAS_SPH)
#ifdef DENSITY_INDEPENDENT_SPH
U = P[i].Mass*SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity / a3, GAMMA_MINUS1);
#else
U = P[i].Mass*SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
else
U = P[i].Mass*SphP[i].Entropy;
#else /* MULTIPHASE */
#ifdef DENSITY_INDEPENDENT_SPH
U = P[i].Mass*SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity / a3, GAMMA_MINUS1);
#else
U = P[i].Mass*SphP[i].Entropy / (GAMMA_MINUS1) * pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
#endif /*MULTIPHASE*/
#endif /*ISOTHERM_EQS*/
for (k=0;k<3;k++)
{
MassCenter[k] = P[i].Mass * P[i].Pos[k];
MV[k] = P[i].Mass * P[i].Vel[k];
MV2[k] = P[i].Mass * P[i].Vel[k] * P[i].Vel[k];
}
dx = Groups[gid].HeadPos[0] - P[i].Pos[0];
dy = Groups[gid].HeadPos[1] - P[i].Pos[1];
dz = Groups[gid].HeadPos[2] - P[i].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
RadiusMax = sqrt(dx * dx + dy * dy + dz * dz);
}
else
{
N = GroupsGet[i].N;
Mass = GroupsGet[i].Mass;
DensityMax = GroupsGet[i].DensityMax;
DensityMin = GroupsGet[i].DensityMin;
PotentialMax = GroupsGet[i].PotentialMax;
W = GroupsGet[i].W;
U = GroupsGet[i].U;
for (k=0;k<3;k++)
{
MassCenter[k] = GroupsGet[i].MassCenter[k];
MV[k] = GroupsGet[i].MV[k];
MV2[k] = GroupsGet[i].MV2[k];
}
RadiusMax = GroupsGet[i].RadiusMax;
}
Groups[gid].N += N;
Groups[gid].Mass += Mass;
Groups[gid].DensityMax = dmax(Groups[gid].DensityMax,DensityMax);
Groups[gid].DensityMin = dmin(Groups[gid].DensityMin,DensityMin);
Groups[gid].PotentialMax = dmax(Groups[gid].PotentialMax,PotentialMax);
Groups[gid].W += W;
Groups[gid].U += U;
for (k=0;k<3;k++)
{
Groups[gid].MassCenter[k] += MassCenter[k];
Groups[gid].MV[k] += MV[k];
Groups[gid].MV2[k] += MV2[k];
}
Groups[gid].RadiusMax = dmax(Groups[gid].RadiusMax,RadiusMax);
}
/*! Do some final opperations on groups properties
*
*/
void fof_final_operations_on_groups_properties()
{
int gid,k;
for (gid=0;gid<Ngroups;gid++)
{
Groups[gid].K =0;
for (k=0;k<3;k++)
{
Groups[gid].MassCenter[k] = Groups[gid].MassCenter[k]/Groups[gid].Mass;
Groups[gid].K += Groups[gid].MV2[k] - (Groups[gid].MV[k]*Groups[gid].MV[k])/Groups[gid].Mass ;
}
/* kinetic energy */
Groups[gid].K *= 0.5;
/* potential energy*/
Groups[gid].W -= (Groups[gid].N*Groups[gid].PotentialMax); /* add offset */
Groups[gid].W *= 0.5;
Groups[gid].T = (Groups[gid].U/Groups[gid].Mass) * GAMMA_MINUS1*All.mumh/All.Boltzmann;
//Groups[gid].Density = Groups[gid].Mass/( 4/3. * PI * pow(Groups[gid].RadiusMax,3) );
Groups[gid].Density = Groups[gid].DensityMax;
/* virial fraction */
Groups[gid].VirialFraction = -2*(Groups[gid].K +Groups[gid].U )/Groups[gid].W;
}
}
/*! Write groups properties
*
*/
void fof_write_groups_properties()
{
FILE *fd;
char buf[500];
int gid;
float vir1,vir2,Wp;
if (Ngroups>0)
{
if (ThisTask==0)
printf("FoF: write groups properties\n");
sprintf(buf, "%s%s.%d", All.OutputDir, "groups", ThisTask);
fd = fopen(buf,"w");
fprintf(fd,"# HeadID N Headx Heady Headz Mass CMx CMy CMz DensityMax DensityMin RadiusMax K U W Wp Vir1 Vir2\n");
for (gid=0;gid<Ngroups;gid++)
if( Groups[gid].Task == ThisTask )
{
//printf("(%02d)--> a true group (head=%d) mass=%g N=%d\n",ThisTask,Groups[gid].HeadID,Groups[i].Mass,Groups[i].N);
Wp = -All.G * Groups[gid].Mass * Groups[gid].Mass / Groups[gid].RadiusMax;
vir1 = -2*(Groups[gid].K +Groups[gid].U )/Groups[gid].W;
vir2 = -2*(Groups[gid].K +Groups[gid].U )/Wp;
fprintf(fd,"%08d %08d %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g \n", Groups[gid].HeadID,Groups[gid].N,Groups[gid].HeadPos[0],Groups[gid].HeadPos[1],Groups[gid].HeadPos[2] , Groups[gid].Mass , Groups[gid].MassCenter[0],Groups[gid].MassCenter[1],Groups[gid].MassCenter[2],Groups[gid].DensityMax,Groups[gid].DensityMin,Groups[gid].RadiusMax,Groups[gid].K,Groups[gid].U,Groups[gid].W,Wp,vir1,vir2);
}
fclose(fd);
}
}
/*! This function computes the gravitational potential for the particles belonging to a group
* due to particles belonging to the same group.
*/
void fof_compute_groups_potential()
{
#ifdef PMGRID
This function do not work with PMGRID on.
#endif
int i;
int n;
#ifndef NOGRAVITY
long long ntot, ntotleft;
int j, k, level, sendTask, recvTask;
int ndone;
int maxfill, ngrp, place, nexport;
int *nsend, *noffset, *nsend_local, *nbuffer, *ndonelist, *numlist;
double fac;
double t0, t1, tstart, tend;
MPI_Status status;
double r2;
int NumPot;
t0 = second();
if(All.ComovingIntegrationOn)
set_softenings();
if(ThisTask == 0)
{
printf("Start computation of potential for particles belonging to a group...\n");
fflush(stdout);
}
/* for the moment, we do not need the tree */
//tstart = second();
//if(TreeReconstructFlag)
// {
// if(ThisTask == 0)
// printf("Tree construction.\n");
//
// force_treebuild(NumPart);
//
// TreeReconstructFlag = 0;
//
// if(ThisTask == 0)
// printf("Tree construction done.\n");
// }
//tend = second();
//All.CPU_TreeConstruction += timediff(tstart, tend);
/* count the number of particles on which we need to compute potential */
for(n = 0, NumPot = 0; n < N_gas; n++)
{
/* keep only if the particle belong to a group */
//SphP[n].FOF_gidx == ...;
if (SphP[n].FOF_gidx > -1)
{
printf("(%d) gidx = %d gid = %d P.ID = %d\n",ThisTask,SphP[n].FOF_gidx,SphP[n].FOF_gid,P[n].ID);
NumPot++;
}
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumPot, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += 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; /* beginn with this index */
ntotleft = ntot; /* particles left for all tasks together */
if (ThisTask==0)
{
printf("FoF: compute potential for %d%09d particles.\n", (int) (ntot / 1000000000), (int) (ntot % 1000000000));
fflush(stdout);
}
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeForceFOF - NTask; i++)
{
if (SphP[i].FOF_gidx > -1) /* the particle belongs to a group or pseudo-group */
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
#ifndef PMGRID
fof_evaluate_groups_potential(i, 0);
#else
!!! need to add a function here
#endif
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
for(k = 0; k < 3; k++)
FOFGravDataGet[nexport].u.Pos[k] = P[i].Pos[k];
#ifdef UNEQUALSOFTENINGS
FOFGravDataGet[nexport].Type = P[i].Type;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
FOFGravDataGet[nexport].Soft = SphP[i].Hsml;
#endif
#endif
FOFGravDataGet[nexport].FOF_gidx = SphP[i].FOF_gidx;
FOFGravDataGet[nexport].FOF_gid = SphP[i].FOF_gid;
GravDataIndexTable[nexport].Task = j;
GravDataIndexTable[nexport].Index = i;
GravDataIndexTable[nexport].SortIndex = nexport;
nexport++;
nsend_local[j]++;
}
}
}
}
qsort(GravDataIndexTable, nexport, sizeof(struct FOFgravdata_index), grav_fof_compare_key);
for(j = 0; j < nexport; j++)
FOFGravDataIn[j] = FOFGravDataGet[FOFGravDataIndexTable[j].SortIndex];
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
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.BunchSizeForceFOF)
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(&FOFGravDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFgravdata_in), MPI_BYTE,
recvTask, TAG_POTENTIAL_A,
&FOFGravDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFgravdata_in), MPI_BYTE,
recvTask, TAG_POTENTIAL_A, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
for(j = 0; j < nbuffer[ThisTask]; j++)
{
#ifndef PMGRID
fof_evaluate_groups_potential(j, 1);
#else
!!! need to add a function here
#endif
}
/* get the result */
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.BunchSizeForceFOF)
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(&FOFGravDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFgravdata_in),
MPI_BYTE, recvTask, TAG_POTENTIAL_B,
&FOFGravDataOut[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFgravdata_in),
MPI_BYTE, recvTask, TAG_POTENTIAL_B, MPI_COMM_WORLD, &status);
/* add the result to the particles */
for(j = 0; j < nsend_local[recvTask]; j++)
{
place = FOFGravDataIndexTable[noffset[recvTask] + j].Index;
SphP[place].FOFPotential += FOFGravDataOut[j + noffset[recvTask]].u.FOFPotential;
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
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);
/* add correction to exclude self-potential */
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
{
/* remove self-potential */
SphP[i].FOFPotential += P[i].Mass / All.SofteningTable[P[i].Type];
if(All.ComovingIntegrationOn)
if(All.PeriodicBoundariesOn)
SphP[i].FOFPotential -= 2.8372975 * pow(P[i].Mass, 2.0 / 3) *
pow(All.Omega0 * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G), 1.0 / 3);
}
/* multiply with the gravitational constant */
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
SphP[i].FOFPotential *= All.G;
//#ifdef PMGRID
//
//need to addapt this...
//
//#ifdef PERIODIC
// pmpotential_periodic();
//#ifdef PLACEHIGHRESREGION
// pmpotential_nonperiodic(1);
//#endif
//#else
// pmpotential_nonperiodic(0);
//#ifdef PLACEHIGHRESREGION
// pmpotential_nonperiodic(1);
//#endif
//#endif
//
//#endif //PMGRID
if(All.ComovingIntegrationOn)
{
#ifndef PERIODIC
fac = -0.5 * All.Omega0 * All.Hubble * All.Hubble;
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
{
for(k = 0, r2 = 0; k < 3; k++)
r2 += P[i].Pos[k] * P[i].Pos[k];
SphP[i].FOFPotential += fac * r2;
}
#endif
}
else
{
#ifndef FORCE_OMEGALAMBDA_TO_0
fac = -0.5 * All.OmegaLambda * All.Hubble * All.Hubble;
if(fac != 0)
{
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
{
for(k = 0, r2 = 0; k < 3; k++)
r2 += P[i].Pos[k] * P[i].Pos[k];
SphP[i].FOFPotential += fac * r2;
}
}
#endif
}
if(ThisTask == 0)
{
printf("potential due to particles in the group done.\n");
fflush(stdout);
}
t1 = second();
All.CPU_Potential += timediff(t0, t1);
#else
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
SphP[i].FOFPotential = 0;
#endif //NOGRAVITY
for(i = 0; i < N_gas; i++)
if (SphP[i].FOF_gidx > -1) /* the particle belong to a group */
{
printf("(%d)--------> pot=%g fofpot=%g (ID=%d)\n",ThisTask,P[i].Potential,SphP[i].FOFPotential,P[i].ID);
}
}
/*! This routine computes the gravitational potential by walking the
* tree. The same opening criteria is used as for the gravitational force
* walk.
*/
void fof_evaluate_groups_potential(int target, int mode)
{
int ptype,next,gidx,gid,i;
double r2, dx, dy, dz, mass, r, u, h, h_inv, wp;
double pot, pos_x, pos_y, pos_z;
#if defined(UNEQUALSOFTENINGS) && !defined(ADAPTIVE_GRAVSOFT_FORGAS)
int maxsofttype;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft = 0;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
pot = 0;
if(mode == 0)
{
pos_x = P[target].Pos[0];
pos_y = P[target].Pos[1];
pos_z = P[target].Pos[2];
ptype = P[target].Type;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = SphP[target].Hsml;
#endif
gidx = SphP[target].FOF_gidx; /* index to the local group */
gid = SphP[target].FOF_gid; /* id of the local group */
}
else
{
pos_x = FOFGravDataGet[target].u.Pos[0];
pos_y = FOFGravDataGet[target].u.Pos[1];
pos_z = FOFGravDataGet[target].u.Pos[2];
#ifdef UNEQUALSOFTENINGS
ptype = FOFGravDataGet[target].Type;
#else
ptype = P[0].Type;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
if(ptype == 0)
soft = FOFGravDataGet[target].Soft;
#endif
gid = FOFGravDataGet[target].FOF_gid;
/* here, gidx = FOFGravDataGet[target].FOF_gidx;
/* do not corresponds to the correct Group index,
/* we need to loop over groups to find the correct one */
for (i=0;i<Ngroups;i++)
{
if ( Groups[i].HeadID == gid) ;
{
gidx = i;
break;
}
if (i==Ngroups-1)
{
printf("(%d) it seems that we have been unable to find the correct group for this Task !",ThisTask);
printf("its better to stop\n");
endrun(1975);
}
}
}
next=Groups[gidx].LocalHead; /* indice towards Head or Pseudo Head */
while(next!=-1) /* loop over all local members */
{
dx = P[next].Pos[0] - pos_x;
dy = P[next].Pos[1] - pos_y;
dz = P[next].Pos[2] - pos_z;
mass = P[next].Mass;
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
h = soft;
if(h < SphP[next].Hsml)
h = SphP[next].Hsml;
#else
h = All.ForceSoftening[ptype]; /* all gas particles are gas particles */
#endif
#endif
r = sqrt(r2);
if(r >= h)
pot -= mass / r;
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
#endif
u = r * h_inv;
if(u < 0.5)
wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
else
wp =
-3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
pot += mass * h_inv * wp;
}
#ifdef PERIODIC
pot += mass * ewald_pot_corr(dx, dy, dz);
#endif
next=SphP[next].FOF_Next;
}
/*
now, we need to fill Exportflag[j], to export the particle if needed.
*/
if(mode == 0)
SphP[target].FOFPotential = pot;
else
FOFGravDataResult[target].u.FOFPotential = pot;
}
/*! Here, each group is flaged to form stars or not.
* The decision is based on physical properties of the group.
*
*/
void fof_set_groups_for_sfr()
{
int gid;
if (ThisTask==0)
printf("FoF: set groups for sfr\n");
Ntsfrgroups=0;
for (gid=0;gid<Ngroups;gid++)
if( Groups[gid].Task == ThisTask ) /* a true group */
{
/* init to 0 */
Groups[gid].SfrFlag = 0;
/* here, we needs to flag for sfr, according to physical quantities */
printf("FoF: (%d) Members=%d MinMembers=%d\n",ThisTask,Groups[gid].N,All.FoF_MinGroupMembers);
if (Groups[gid].N<All.FoF_MinGroupMembers)
continue;
printf("FoF: (%d) Mass=%g MinMass%g\n",ThisTask,Groups[gid].Mass*All.CMUtoMsol,All.FoF_MinGroupMass*All.CMUtoMsol);
if (Groups[gid].Mass<All.FoF_MinGroupMass)
continue;
#ifndef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
/* remove group that are unbound 2*(T+U)/W >~ 1 */
printf("FoF: (%d) VirialFraction = %g FoF_VirialFractionThreshold = %g\n",ThisTask,Groups[gid].VirialFraction,All.FoF_VirialFractionThreshold);
if (Groups[gid].VirialFraction > All.FoF_VirialFractionThreshold)
continue;
#endif
#ifdef SFR
printf("FoF: (%d) Temperature=%g MaxTemperature%g\n",ThisTask, Groups[gid].T, All.StarFormationTemperature);
if (Groups[gid].T>All.StarFormationTemperature)
continue;
/* compute the sfr probability */
double p,r;
double dt=0;
double tstar=1.;
double MeanBaryonicDensity;
switch (All.StarFormationType)
{
case 0:
All.ThresholdDensity = All.StarFormationDensity;
break;
case 1:
MeanBaryonicDensity = All.OmegaBaryon * (3 * All.Hubble * All.Hubble / (8 * M_PI * All.G));
All.ThresholdDensity = dmax(All.StarFormationDensity/a3inv,200*MeanBaryonicDensity);
break;
case 2:
All.ThresholdDensity = All.StarFormationDensity;
All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
break;
case 3:
All.StarFormationTime = 1./All.StarFormationCstar / pow(4*PI*All.G*All.StarFormationDensity,0.5);
/* All.ThresholdDensity is fixed further */
break;
case 4:
All.ThresholdDensity = All.StarFormationDensity;
All.StarFormationTime = 1./All.StarFormationCstar / pow((32/(3*PI))*All.G*All.StarFormationDensity,0.5);
break;
case 6:
/* estimate the freefall time */
All.StarFormationTime = 1./ pow((32/(3*PI))*All.G*Groups[gid].Density*a3inv,0.5);
break;
}
dt = DeltaTime; /* computed earlier */
//tstar = All.StarFormationTime / pow(Groups[gid].Density*a3inv/All.StarFormationDensity,0.5);
tstar = All.StarFormationTime * All.StarFormationCstar;
p = (1. - exp(-dt/tstar) );
if (All.FoF_StarFormationType==1)
p = 1;
r = get_StarFormation_random_number(Groups[gid].HeadID);
printf("FoF: (%d) Probability=%g Random=%g tstar=%g dt=%g GroupDensity=%g StarFormationDensity=%g \n",ThisTask, p,r,tstar,dt,Groups[gid].Density,All.StarFormationDensity );
if (r>p)
continue;
#endif
/* the groups may form stars */
Groups[gid].SfrFlag = 1;
Ntsfrgroups+=1;
//printf("%08d %08d %g %g %g %g %g %g %g %g %g %g %g %g %g \n", Groups[gid].HeadID,Groups[gid].N,Groups[gid].HeadPos[0],Groups[gid].HeadPos[1],Groups[gid].HeadPos[2] , Groups[gid].Mass , Groups[gid].MassCenter[0],Groups[gid].MassCenter[1],Groups[gid].MassCenter[2],Groups[gid].DensityMax,Groups[gid].DensityMin,Groups[gid].RadiusMax,Groups[gid].K,Groups[gid].U,Groups[gid].W);
}
MPI_Allreduce(&Ntsfrgroups, &Tot_Ntsfrgroups, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (ThisTask==0)
printf("FoF: %d true groups flagged for star formation.\n",Tot_Ntsfrgroups);
}
#ifdef SFR
/*! Flag particles for star formation based on their groups
*
*/
void fof_star_formation()
{
int gid;
int next;
int i;
FILE *fd;
char buf[500];
double mnewstars=0,Tot_mnewstars=0;
double mstars=0,Tot_mstars=0;
int Nnewstars=0,Tot_Nnewstars=0;
double star_formation_rate;
if (ThisTask==0)
printf("FoF: star formation\n");
if (Ngroups>0)
{
sprintf(buf, "%s%s.%d", All.OutputDir, "sfp", ThisTask);
fd = fopen(buf,"w");
for (gid=0;gid<Ngroups;gid++)
{
if( Groups[gid].SfrFlag == 1 ) /* the group must be flagged as forming stars */
{
/* loop over particles of the group */
next= Groups[gid].LocalHead ;
while(next!=-1)
{
/* create a new star */
P[next].Type = ST;
RearrangeParticlesFlag=1;
/* count mass */
mnewstars+=P[i].Mass;
Nnewstars+=1;
fprintf(fd,"%d\n",P[next].ID);
next=SphP[next].FOF_Next;
}
}
}
fflush(fd);
fclose(fd);
}
/* compute star mass */
mstars=0;
for (i=0;i< NumPart;i++)
{
if (P[i].Type==ST)
mstars+= P[i].Mass;
}
/* share results */
MPI_Allreduce(&mstars, &Tot_mstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&mnewstars, &Tot_mnewstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&Nnewstars, &Tot_Nnewstars, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (All.TimeStep>0)
star_formation_rate = Tot_mnewstars/All.TimeStep;
else
star_formation_rate = 0.;
#ifndef PY_INTERFACE
if (ThisTask==0)
{
//fprintf(FdSfr, "Step %d, Time: %g, NewStarsPart: %g \n", All.NumCurrentTiStep, All.Time, Tot_mstars);
fprintf(FdSfr, "%15g %15g %15g %15g %15g\n", All.Time,All.TimeStep,Tot_mstars,Tot_mnewstars,star_formation_rate);
fflush(FdSfr);
}
#endif
if(ThisTask == 0)
{
printf("FoF: Total number of new star particles : %d \n",(int) Tot_Nnewstars);
printf("FoF: Total number of star particles : %d%09d\n",(int) (All.TotN_stars / 1000000000), (int) (All.TotN_stars % 1000000000));
fflush(stdout);
}
#ifdef CHECK_BLOCK_ORDER
int typenow;
rearrange_particle_sequence();
typenow = 0;
for(i = 0; i < NumPart; i++)
{
if ((P[i].Type<typenow)&&(typenow<=ST))
{
printf("\nBlock order error\n");
printf("(%d) i=%d id=%d Type=%d typenow=%d\n\n",ThisTask,i,P[i].ID,P[i].Type,typenow);
endrun(999004);
}
typenow = P[i].Type;
}
#endif
#ifdef CHECK_ID_CORRESPONDENCE
rearrange_particle_sequence();
for(i = N_gas; i < N_gas+N_stars; i++)
{
if( StP[P[i].StPIdx].PIdx != i )
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in starformation) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
endrun(999005);
}
if(StP[P[i].StPIdx].ID != P[i].ID)
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in starformation) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
endrun(999006);
}
}
#endif
/* force domaine decomposition */
if (Tot_Nnewstars>1)
All.NumForcesSinceLastDomainDecomp = All.TotNumPart * All.TreeDomainUpdateFrequency *10;
}
#endif /* SFR */
/*! This function is responsible of sampling
* the IMF.
*/
void fof_IMF_sampling(int Np)
{
int i,j,k;
double fp;
char buf[500];
FILE *fd = 0;
double Mcl=0;
for (i=0;i<Np;i++)
{
Mcl += FoF_P[i].Mass*All.CMUtoMsol;
}
if (ThisTask==0)
printf("\nFoF: fof_IMF_sampling : members = %6d cluster mass = %6g [Msol]\n",Np,Mcl);
#ifdef CHIMIE_OPTIMAL_SAMPLING
double current_mass;
double Mtot;
int stop_discrete_imf;
double Mrem_cl;
double Mrem_p;
double m1,m2,mp;
current_mass = optimal_init_norm(Mcl); /* in Msol, compute Cp->k and Cp->m_max */
/*
* do the descrete part of the IMF
*/
stop_discrete_imf=0;
for (k=0;k<Np;k++)
{
if (stop_discrete_imf)
break;
/* add at least one star */
Mtot+=current_mass;
FoF_P[k].MassNew = current_mass;
FoF_P[k].MassMax = current_mass;
FoF_P[k].MassMin = current_mass;
FoF_P[k].MassSSP = Mcl;
FoF_P[k].NStars = 1;
FoF_P[k].k = optimal_get_k_normalisation_factor(); /* imf normalisation, Cp->k is set by optimal_init_norm above */
if (optimal_stop_loop(current_mass))
break;
current_mass = optimal_get_next_mass(current_mass);
while( fabs( FoF_P[k].MassNew-FoF_P[k].Mass*All.CMUtoMsol + current_mass ) < fabs( FoF_P[k].MassNew-FoF_P[k].Mass*All.CMUtoMsol ) )
{
Mtot+=current_mass;
FoF_P[k].MassNew += current_mass;
FoF_P[k].MassMin = current_mass;
FoF_P[k].NStars++;
if ( FoF_P[k].NStars >= MAX_STARS_PER_PART)
{
stop_discrete_imf=1;
break;
}
if (optimal_stop_loop(current_mass))
break;
current_mass = optimal_get_next_mass(current_mass);
}
}
/*
* do the continuous part of the IMF
*/
// remaining mass in the cluster
Mrem_cl = Mcl - Mtot;
// mass of remaining particles
Mrem_p=0;
j = k;
for (k=j;k<Np;k++)
Mrem_p+=FoF_P[k].Mass*All.CMUtoMsol;
for (k=j;k<Np;k++)
{
mp = (FoF_P[k].Mass*All.CMUtoMsol)*Mrem_cl/Mrem_p;
m2 = FoF_P[k-1].MassMin;
m1 = optimal_get_m1_from_m2(m2,mp);
Mtot += mp;
FoF_P[k].MassNew = mp;
FoF_P[k].MassMax = m2;
FoF_P[k].MassMin = m1;
FoF_P[k].NStars = -1;
FoF_P[k].MassSSP = Mcl;
}
/* write some info */
if (ThisTask==0)
{
sprintf(buf, "%s%s_%04d.txt", All.OutputDir, All.FoF_IMFSnapshotFileBase, All.FoF_IMFSnapshotFileCount++);
if(!(fd = fopen(buf, "w")))
{
printf("can't open file `%s' for savinf fof.\n", buf);
}
else
{
fprintf(fd,"fof_IMF_sampling :Time=%g members = %6d cluster mass = %6g [Msol]\n",All.Time,Np,Mcl);
fprintf(fd,"\n i ID #stars MassMin MassMax CurrentMass MassNew MassRatio\n\n");
for (i=0;i<Np;i++)
{
fp = 100*fabs(FoF_P[i].Mass*All.CMUtoMsol - FoF_P[i].MassNew) / (FoF_P[i].Mass*All.CMUtoMsol);
//if (fp>50)
//printf("FoF: fof_IMF_sampling : %3d %6d %4d %8.4g %8.4g %6.2f %6.2f %g\n",i,FoF_P[i].ID,FoF_P[i].NStars,FoF_P[i].MassMin,FoF_P[i].MassMax,FoF_P[i].MassNew,FoF_P[i].Mass*All.CMUtoMsol, FoF_P[i].Mass*All.CMUtoMsol/FoF_P[i].MassNew );
fprintf(fd,"%3d %6d %4d %8.4g %8.4g %6.2f %6.2f %6.3f %g %g %g\n",i, FoF_P[i].ID, FoF_P[i].NStars, FoF_P[i].MassMin, FoF_P[i].MassMax, FoF_P[i].Mass*All.CMUtoMsol, FoF_P[i].MassNew, FoF_P[i].Mass*All.CMUtoMsol/FoF_P[i].MassNew, FoF_P[i].Pos[0], FoF_P[i].Pos[1], FoF_P[i].Pos[2] );
}
/* check */
fp = 100* (Mcl-Mtot)/(Mcl);
//if (fp>1e-3)
{
printf("FoF: fof_IMF_sampling : Mtot=%8.2f mass difference %7.3f %\n", Mtot,100* (Mcl-Mtot)/(Mcl) );
fprintf(fd,"\nMtot=%8.2f mass difference %7.3f %\n", Mtot,100* (Mcl-Mtot)/(Mcl) );
}
fflush(fd);
}
fclose(fd);
}
/* convert all masses in Code Mass Unit */
for (k=0;k<Np;k++)
{
FoF_P[k].MassMax *= All.MsoltoCMU;
FoF_P[k].MassMin *= All.MsoltoCMU;
FoF_P[k].MassSSP *= All.MsoltoCMU;
}
#endif
}
#ifdef SFR
/*! Flag particles for star formation based on their groups
* in addition, we address to each particles Min and Max stellar
* masses that it will represent, according to an IMF
*
*/
void fof_star_formation_and_IMF_sampling()
{
int gid;
int next;
int index;
int i,j,k;
double mnewstars=0,Tot_mnewstars=0;
double mstars=0,Tot_mstars=0;
int Nnewstars=0,Tot_Nnewstars=0;
double star_formation_rate;
int ntop;
int *ntoplist, *ntopoffset;
int *npartlist, *npartoffset;
int Np=0,Np_local;
int Ns=0;
int pass;
if (ThisTask==0)
printf("FoF: star formation and IMF sampling\n");
/* create the TopTGroups structure
TopTSfrGroups is a structure shared by all Tasks and contains
the Task and Head of all true groups.
*/
TopTSfrGroups = malloc(sizeof(struct fof_topgroups_exchange)*Tot_Ntsfrgroups);
TopTSfrGroups_local = malloc(sizeof(struct fof_topgroups_exchange)*Ntsfrgroups);
j=0;
if (Ngroups>0)
{
for (gid=0;gid<Ngroups;gid++)
{
if( Groups[gid].Task == ThisTask ) /* a true group */
{
if( Groups[gid].SfrFlag == 1 ) /* the group must be flagged as forming stars */
{
if (j==Ntsfrgroups)
{
printf("j=%d == Nsfrgroups=%d !!!\n",j,Ntsfrgroups);
endrun(1970);
}
TopTSfrGroups_local[j].Task = Groups[gid].Task;
TopTSfrGroups_local[j].Head = Groups[gid].Head;
j++;
}
}
}
}
ntoplist = malloc(sizeof(int) * NTask);
ntopoffset = malloc(sizeof(int) * NTask);
MPI_Allgather(&Ntsfrgroups, 1, MPI_INT, ntoplist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntopoffset[0] = 0; i < NTask; i++)
{
if(i > 0)
ntopoffset[i] = ntopoffset[i - 1] + ntoplist[i - 1];
}
for(i = 0; i < NTask; i++)
{
ntoplist[i] *= sizeof(struct fof_topgroups_exchange);
ntopoffset[i] *= sizeof(struct fof_topgroups_exchange);
}
MPI_Allgatherv(TopTSfrGroups_local, Ntsfrgroups * sizeof(struct fof_topgroups_exchange), MPI_BYTE, TopTSfrGroups, ntoplist, ntopoffset, MPI_BYTE, MPI_COMM_WORLD);
free(ntoplist);
free(ntopoffset);
/* check that TopTGroups is correct */
/*
for (i=0;i<Ntsfrgroups;i++)
printf(">>> local (%d) i=%d Task=%d Head=%d \n",ThisTask,i,TopTSfrGroups_local[i].Task,TopTSfrGroups_local[i].Head);
if (ThisTask==0)
for (i=0;i<Tot_Ntsfrgroups;i++)
printf(">>> global (%d) i=%d Task=%d Head=%d \n",ThisTask,i,TopTSfrGroups[i].Task,TopTSfrGroups[i].Head);
*/
/*
Here, we want to fill the FoF_P structure that contains particles information for a group that will form an IMF. This structure is shared by all Tasks.
Loop over all TopTSfrGroups and check if a local Group corresponds to it.
1) fill FoF_P_local with local information
2) send FoF_P_local to all and create FoF_P
*/
RearrangeParticlesFlag=0;
for (i=0;i<Tot_Ntsfrgroups;i++)
{
Np_local=-1;
pass=0;
if (Ngroups>0)
{
for (gid=0;gid<Ngroups;gid++)
{
if( Groups[gid].SfrFlag == 1 ) /* the group must be flagged as forming stars */
{
if ( (Groups[gid].Task==TopTSfrGroups[i].Task) && ((Groups[gid].Head==TopTSfrGroups[i].Head))) /* the pseudo group fit the top group */
{
/* count number of stars */
/* Note : Groups[gid].N != Np_local as it corresponds to the total group for the true Groups */
Np_local = 0;
next= Groups[gid].LocalHead ;
while(next!=-1)
{
Np_local+=1;
next=SphP[next].FOF_Next;
}
if (pass!=0)
{
printf("(%d) we already found a group linked to Task=%d Head=%d in this proc. \n",ThisTask,TopTSfrGroups[i].Task,TopTSfrGroups[i].Head);
printf("we should reduce such groups before... we prefere to stop.\n");
endrun(1974);
}
pass++;
/* prepare stars */
FoF_P_local = malloc(Np_local*sizeof(struct fof_particles_in));
next= Groups[gid].LocalHead ;
j=0;
while(next!=-1)
{
FoF_P_local[j].Mass = P[next].Mass;
FoF_P_local[j].Task = ThisTask;
FoF_P_local[j].Index = next;
FoF_P_local[j].ID = P[next].ID;
FoF_P_local[j].Pos[0] = P[next].Pos[0];
FoF_P_local[j].Pos[1] = P[next].Pos[1];
FoF_P_local[j].Pos[2] = P[next].Pos[2];
#ifdef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
FoF_P_local[j].Vel[0] = P[next].Vel[0];
FoF_P_local[j].Vel[1] = P[next].Vel[1];
FoF_P_local[j].Vel[2] = P[next].Vel[2];
#ifdef ISOTHERM_EQS
FoF_P_local[j].EnergySpec = SphP[next].Entropy;
#else
#ifdef DENSITY_INDEPENDENT_SPH
FoF_P_local[j].EnergySpec = SphP[next].Entropy / (GAMMA_MINUS1) * pow(SphP[next].EgyWtDensity / a3, GAMMA_MINUS1);
#else
FoF_P_local[j].EnergySpec = SphP[next].Entropy / (GAMMA_MINUS1) * pow(SphP[next].Density / a3, GAMMA_MINUS1);
#endif
#endif /*ISOTHERM_EQS*/
#ifdef DENSITY_INDEPENDENT_SPH
FoF_P_local[j].Density = SphP[next].EgyWtDensity;
#else
FoF_P_local[j].Density = SphP[next].Density;
#endif
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
FoF_P_local[j].Hsml = P[next].Hsml;
#endif
#endif
j++;
next=SphP[next].FOF_Next;
}
/* sort particles according to some physical quantity */
qsort(FoF_P_local, Np_local, sizeof(struct fof_particles_in), fof_compare_particles);
}
}
}
}
/* all Tasks get all particles and sample the IMF, thus, we no longer needs to communicate back */
if (Np_local==-1) /* no particles for this Task */
{
Np_local = 0;
FoF_P_local = malloc(Np_local*sizeof(struct fof_particles_in));
}
npartlist = malloc(sizeof(int) * NTask);
npartoffset = malloc(sizeof(int) * NTask);
MPI_Allgather(&Np_local, 1, MPI_INT, npartlist, 1, MPI_INT, MPI_COMM_WORLD);
for(k = 0, Np=0, npartoffset[0] = 0; k < NTask; k++)
{
Np += npartlist[k];
if(k > 0)
npartoffset[k] = npartoffset[k - 1] + npartlist[k - 1];
}
FoF_P = malloc(Np * sizeof(struct fof_particles_in));
for(k = 0; k < NTask; k++)
{
npartlist[k] *= sizeof(struct fof_particles_in);
npartoffset[k] *= sizeof(struct fof_particles_in);
}
MPI_Allgatherv(FoF_P_local, Np_local * sizeof(struct fof_particles_in), MPI_BYTE, FoF_P, npartlist, npartoffset, MPI_BYTE, MPI_COMM_WORLD);
/* sort particles according to some physical quantity */
qsort(FoF_P, Np, sizeof(struct fof_particles_in), fof_compare_particles);
//Ns += Np;
//if (ThisTask==0)
// printf("check i=%d Np=%d Ns=%d Task=%d Head=%d\n",i,Np,Ns,TopTSfrGroups[i].Task,TopTSfrGroups[i].Head);
//if (ThisTask==0)
//{
// printf("check %d \n",i);
// for(k=0;k<Np;k++)
// printf("check %d %g\n",FoF_P[k].ID,FoF_P[k].Mass);
//}
#ifdef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
double VirialFraction;
double Epot,Ekin,Eint;
VirialFraction = fof_compute_groups_virial_locally(Np,&Epot,&Ekin,&Eint);
if (ThisTask==0)
fof_save_fof_particles(Np,VirialFraction,Epot,Ekin,Eint);
if (VirialFraction < All.FoF_VirialFractionThreshold)
{
#endif
/*
now, sample the IMF
*/
fof_IMF_sampling(Np);
/* now, each proc copy its particles from FoF_P */
for(k = 0; k < Np; k++)
{
if (FoF_P[k].Task==ThisTask)
{
index = FoF_P[k].Index;
/* at this point, we turn gas particles into stellar particles */
#ifdef STELLAR_PROP
if (N_stars+Nnewstars+1 > All.MaxPartStars)
{
printf("No memory left for new star particle : N_stars+1=%d MaxPartStars=%d\n\n",N_stars+1,All.MaxPartStars);
endrun(999008);
}
#else
if (N_stars+Nnewstars+1 > All.MaxPart-N_gas)
{
printf("No memory left for new star particle : N_stars+1=%d All.MaxPart-N_gas=%d\n\n",N_stars+1,All.MaxPart-N_gas);
endrun(999009);
}
#endif
RearrangeParticlesFlag=1;
P[index].Mass = FoF_P[k].MassNew*All.MsoltoCMU;
P[index].Type = ST;
//P[index].FoFidx = k; /* reference to the FoF_P structure (no longer used) */
/* we store IMF information in the SPH structure in order to use the information in the rearrange_particle_sequence routine
and as StP still doesn't exists for these particles. */
SphP[index].FOF_MassMax = FoF_P[k].MassMax;
SphP[index].FOF_MassMin = FoF_P[k].MassMin;
SphP[index].FOF_MassSSP = FoF_P[k].MassSSP;
SphP[index].FOF_NStars = FoF_P[k].NStars;
#ifdef CHIMIE_OPTIMAL_SAMPLING
SphP[index].OptIMF_k = FoF_P[k].k;
SphP[index].OptIMF_CurrentMass = FoF_P[k].MassMax*All.CMUtoMsol;
#endif
mnewstars+=P[index].Mass;
Nnewstars+=1;
/* gather energy */
#ifdef FEEDBACK
LocalSysState.StarEnergyFeedback += SphP[i].EgySpecFeedback * P[i].Mass;
#endif
}
}
#ifdef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
}
#endif
free(npartlist);
free(npartoffset);
free(FoF_P);
free(FoF_P_local);
}
free(TopTSfrGroups_local);
free(TopTSfrGroups);
/* rearrange_particle_sequence
* this may not be done after each group, as we
* need to conserve particle order
*/
rearrange_particle_sequence();
/* compute star mass */
mstars=0;
for (i=0;i< NumPart;i++)
{
if (P[i].Type==ST)
mstars+= P[i].Mass;
}
/* share results */
MPI_Allreduce(&mstars, &Tot_mstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&mnewstars, &Tot_mnewstars, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(&Nnewstars, &Tot_Nnewstars, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
if (All.TimeStep>0)
star_formation_rate = Tot_mnewstars/All.TimeStep;
else
star_formation_rate = 0.;
#ifndef PY_INTERFACE
if (ThisTask==0)
{
//fprintf(FdSfr, "Step %d, Time: %g, NewStarsPart: %g \n", All.NumCurrentTiStep, All.Time, Tot_mstars);
fprintf(FdSfr, "%15g %15g %15g %15g %15g\n", All.Time,All.TimeStep,Tot_mstars,Tot_mnewstars,star_formation_rate);
fflush(FdSfr);
}
#endif
if(ThisTask == 0)
{
printf("FoF: Total number of new star particles : %d \n",(int) Tot_Nnewstars);
printf("FoF: Total number of star particles : %d%09d\n",(int) (All.TotN_stars / 1000000000), (int) (All.TotN_stars % 1000000000));
fflush(stdout);
}
#ifdef CHECK_BLOCK_ORDER
int typenow;
rearrange_particle_sequence();
typenow = 0;
for(i = 0; i < NumPart; i++)
{
if ((P[i].Type<typenow)&&(typenow<=ST))
{
printf("\nBlock order error\n");
printf("(%d) i=%d id=%d Type=%d typenow=%d\n\n",ThisTask,i,P[i].ID,P[i].Type,typenow);
endrun(999004);
}
typenow = P[i].Type;
}
#endif
#ifdef CHECK_ID_CORRESPONDENCE
rearrange_particle_sequence();
for(i = N_gas; i < N_gas+N_stars; i++)
{
if( StP[P[i].StPIdx].PIdx != i )
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in starformation) N_stars=%d N_gas=%d i=%d id=%d P[i].StPIdx=%d StP[P[i].StPIdx].PIdx=%d\n\n",ThisTask,N_stars,N_gas,i,P[i].ID,P[i].StPIdx,StP[P[i].StPIdx].PIdx);
endrun(999005);
}
if(StP[P[i].StPIdx].ID != P[i].ID)
{
printf("\nP/StP correspondance error\n");
printf("(%d) (in starformation) N_gas=%d N_stars=%d i=%d Type=%d P.Id=%d P[i].StPIdx=%d StP[P[i].StPIdx].ID=%d \n\n",ThisTask,N_gas,N_stars,i,P[i].Type,P[i].ID, P[i].StPIdx, StP[P[i].StPIdx].ID);
endrun(999006);
}
}
#endif
/* force domaine decomposition */
if (Tot_Nnewstars>1)
All.NumForcesSinceLastDomainDecomp = All.TotNumPart * All.TreeDomainUpdateFrequency *10;
}
#endif /* SFR */
#ifdef FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
/*! This function is responsible of sampling
* the IMF.
*/
double fof_compute_groups_virial_locally(int Np,double *Epot,double *Ekin,double *Eint)
{
int i,j,k;
double U,K,W;
double pot;
double r2, dx, dy, dz, massi,massj, r, u, h, h_inv, wp;
double pos_x,pos_y,pos_z;
double fac;
double VirialFraction;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
double soft;
#endif
#ifdef PERIODIC
double boxsize, boxhalf;
boxsize = All.BoxSize;
boxhalf = 0.5 * All.BoxSize;
#endif
if (ThisTask==0)
printf("\n(%d)FoF: fof_compute_groups_virial_locally : members = %6d \n",ThisTask,Np);
/* compute potential energy */
for (i=0;i<Np;i++)
FoF_P[i].Potential=0;
for (i=0;i<Np;i++)
{
pos_x = FoF_P[i].Pos[0];
pos_y = FoF_P[i].Pos[1];
pos_z = FoF_P[i].Pos[2];
massi = FoF_P[i].Mass;
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
soft = FoF_P[i].Hsml;
#endif
for (j=i+1;j<Np;j++)
{
dx = FoF_P[j].Pos[0] - pos_x;
dy = FoF_P[j].Pos[1] - pos_y;
dz = FoF_P[j].Pos[2] - pos_z;
massj = FoF_P[j].Mass;
#ifdef PERIODIC
dx = NEAREST(dx);
dy = NEAREST(dy);
dz = NEAREST(dz);
#endif
r2 = dx * dx + dy * dy + dz * dz;
#ifdef UNEQUALSOFTENINGS
#ifdef ADAPTIVE_GRAVSOFT_FORGAS
h = soft;
if(h < FoF_P[j].Hsml)
h = FoF_P[j].Hsml;
#else
h = All.ForceSoftening[0]; /* all gas particles are gas particles */
#endif
#endif
r = sqrt(r2);
if(r >= h)
pot -= 1 / r;
else
{
#ifdef UNEQUALSOFTENINGS
h_inv = 1.0 / h;
#endif
u = r * h_inv;
if(u < 0.5)
wp = -2.8 + u * u * (5.333333333333 + u * u * (6.4 * u - 9.6));
else
wp = -3.2 + 0.066666666667 / u + u * u * (10.666666666667 +
u * (-16.0 + u * (9.6 - 2.133333333333 * u)));
pot = h_inv * wp;
}
#ifdef PERIODIC
pot += ewald_pot_corr(dx, dy, dz);
#endif
FoF_P[i].Potential+=pot*massj;
FoF_P[j].Potential+=pot*massi;
}
}
/* add some corrections */
for (i=0;i<Np;i++)
{
if(All.ComovingIntegrationOn)
if(All.PeriodicBoundariesOn)
FoF_P[i].Potential -= 2.8372975 * pow(FoF_P[i].Mass, 2.0 / 3) *
pow(All.Omega0 * 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G), 1.0 / 3);
FoF_P[i].Potential *= All.G;
if(All.ComovingIntegrationOn)
{
#ifndef PERIODIC
fac = -0.5 * All.Omega0 * All.Hubble * All.Hubble;
for(i = 0; i < N_gas; i++)
{
for(k = 0, r2 = 0; k < 3; k++)
r2 += FoF_P[i].Pos[k] * FoF_P[i].Pos[k];
FoF_P[i].Potential += fac * r2;
}
#endif
}
else
{
#ifndef FORCE_OMEGALAMBDA_TO_0
fac = -0.5 * All.OmegaLambda * All.Hubble * All.Hubble;
if(fac != 0)
{
for(i = 0; i < N_gas; i++)
{
for(k = 0, r2 = 0; k < 3; k++)
r2 += FoF_P[i].Pos[k] * FoF_P[i].Pos[k];
FoF_P[i].Potential += fac * r2;
}
}
#endif
}
}
W = 0;
for (i=0;i<Np;i++)
W+= FoF_P[i].Mass * FoF_P[i].Potential;
W = 0.5*W;
/* compute kinetic energy */
K = 0;
for (i=0;i<Np;i++)
K+= FoF_P[i].Mass * ( FoF_P[i].Vel[0]*FoF_P[i].Vel[0] + FoF_P[i].Vel[1]*FoF_P[i].Vel[1] + FoF_P[i].Vel[2]*FoF_P[i].Vel[2]);
K = 0.5*K;
/* compute internal energy */
U = 0;
for (i=0;i<Np;i++)
{
if (ThisTask==0)
printf("--> %d %g %g\n",FoF_P[i].ID,FoF_P[i].Mass,FoF_P[i].EnergySpec);
U+= FoF_P[i].Mass * FoF_P[i].EnergySpec;
}
VirialFraction=-2*(K + U )/ W;
if (ThisTask==0)
printf("\n(%d)FoF: Time=%15g VirialFraction = %15g \n",ThisTask,All.Time,VirialFraction);
*Epot = W;
*Ekin = K;
*Eint = U;
return VirialFraction;
}
void fof_save_fof_particles(int Np,double VirialFraction,double Epot,double Ekin,double Eint)
{
FILE *fd;
char buf[500];
int i;
if (VirialFraction < All.FoF_VirialFractionThreshold)
sprintf(buf, "%s%s_accepted_%04d.txt", All.OutputDir, All.FoF_SnapshotFileBase, All.FoF_SnapshotFileCount_Accepted++);
else
//sprintf(buf, "%s%s_rejected_%04d.txt", All.OutputDir, All.FoF_SnapshotFileBase, All.FoF_SnapshotFileCount_Rejected++);
return;
printf("FoF: fof_save_fof_particles\n");
fd = fopen(buf, "w");
fprintf(fd,"Time = %g\nMembers = %6d\nVirialFraction = %g\nEpot = %g\nEkin = %g\nEint = %g\n\n",All.Time,Np,VirialFraction,Epot,Ekin,Eint);
for (i=0;i<Np;i++)
{
fprintf(fd,"%3d %20.10e %20.15e %20.10e %20.10e %20.10e %20.10e %20.10e %20.10e %20.10e\n",FoF_P[i].ID,FoF_P[i].Pos[0],FoF_P[i].Pos[1],FoF_P[i].Pos[2],FoF_P[i].Vel[0],FoF_P[i].Vel[1],FoF_P[i].Vel[2],FoF_P[i].Mass,FoF_P[i].EnergySpec,FoF_P[i].Density);
}
fclose(fd);
}
#endif // FOF_COMPUTE_GROUPS_VIRIAL_LOCALLY
void fof_group_info(void)
{
int i;
int next;
for(i = 0; i < N_gas; i++)
if ( (P[i].Type == 0) && (SphP[i].Density*a3inv>All.FoF_ThresholdDensity) )
if (P[i].ID==3480)
if (SphP[i].FOF_Head==i)
{
printf("(%d) this is indeed a head (CPUHead=%d CPUPrev=%d prev=%d)\n",ThisTask,SphP[i].FOF_CPUHead,SphP[i].FOF_CPUPrev,SphP[i].FOF_Prev);
next=i;
while(next!=-1)
{
printf("%08d\n",P[next].ID);
next=SphP[next].FOF_Next;
}
}
}
/*! Send pseudo heads and find correct head in distant CPU
*
* Here, we call the routinie : fof_link_pseudo_heads(j)
*
*/
void fof_send_and_link_pseudo_heads(void)
{
long long ntot, ntotleft;
int i, j, k, n, ngrp, maxfill, source, ndone;
int level, sendTask, recvTask, nexport, place;
int *nbuffer, *noffset, *nsend_local, *nsend, *numlist, *ndonelist;
int npleft;
int iter=0;
MPI_Status status;
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NpHead, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += 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);
do
{
i=0;
ntotleft = ntot; /* particles left for all tasks together */
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeFOF - NTask; i++)
if ( (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead!=ThisTask) && (SphP[i].FOF_Done==0) ) /* a pseudo head */
{
ndone++;
j = SphP[i].FOF_CPUHead;
for(k = 0; k < 3; k++)
{
FOFDataIn[nexport].Pos[k] = P[i].Pos[k];
}
FOFDataIn[nexport].Density = SphP[i].Density;
FOFDataIn[nexport].Hsml = SphP[i].Hsml;
FOFDataIn[nexport].ID = P[i].ID;
FOFDataIn[nexport].FOF_Prev = SphP[i].FOF_Prev;
FOFDataIn[nexport].FOF_Head = SphP[i].FOF_Head;
FOFDataIn[nexport].Index = i;
FOFDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
qsort(FOFDataIn, nexport, sizeof(struct FOFdata_in), fof_compare_key);
for(j = 1, noffset[0] = 0; j < NTask; j++)
noffset[j] = noffset[j - 1] + nsend_local[j - 1];
MPI_Allgather(nsend_local, NTask, MPI_INT, nsend, NTask, MPI_INT, MPI_COMM_WORLD);
/* now do the particles that need to be exported */
for(level = 1; level < (1 << PTask); level++)
{
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.BunchSizeFOF)
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(&FOFDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFdata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&FOFDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFdata_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];
}
/* now do the imported particles */
for(j = 0; j < nbuffer[ThisTask]; j++)
fof_link_pseudo_heads(j);
/* do a block to measure imbalance */
/* get the result */
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.BunchSizeFOF)
break;
sendTask = ThisTask;
recvTask = ThisTask ^ ngrp;
if(recvTask < NTask)
{
if(nsend[ThisTask * NTask + recvTask] > 0 || nsend[recvTask * NTask + ThisTask] > 0)
{
MPI_Sendrecv(&FOFDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct FOFdata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B,
&FOFDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct FOFdata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B, MPI_COMM_WORLD, &status);
for(j = 0; j < nsend_local[recvTask]; j++)
{
source = j + noffset[recvTask];
place = FOFDataIn[source].Index;
SphP[place].FOF_Head = FOFDataPartialResult[source].FOF_Head;
SphP[place].FOF_CPUHead = FOFDataPartialResult[source].FOF_CPUHead;
SphP[place].FOF_Prev = FOFDataPartialResult[source].FOF_Prev;
SphP[place].FOF_CPUPrev = FOFDataPartialResult[source].FOF_CPUPrev;
SphP[place].FOF_Done = FOFDataPartialResult[source].FOF_Done;
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
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];
}
/* do final operations on results */
/* count particles that needs another loop */
MPI_Barrier(MPI_COMM_WORLD);
for(i = 0, npleft = 0; i < N_gas; i++)
if ( (SphP[i].FOF_Head==i) && (SphP[i].FOF_CPUHead!=ThisTask) ) /* a pseudo head */
if (SphP[i].FOF_Done==0) /* not done yet */
npleft++;
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&npleft, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
if(ntot > 0)
{
iter++;
if (ThisTask==0)
{
printf("FoF: iteration %03d : need to repeat for %d%09d particles.\n",iter,(int) (ntot / 1000000000), (int) (ntot % 1000000000));
fflush(stdout);
}
if (iter>10)
{
endrun(1260);
}
}
}
while(ntot>0);
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
}
void fof_allocate_groups(void)
{
Groups = malloc(sizeof(struct fof_groups_data) * (NpHead+NHead));
}
void fof_free_groups(void)
{
free(Groups);
}
/*! Link an imported pseudo head to a local head
*
* If the imported particle points towards a particle
* that is actually exported, wait for another loop..
* and flag it as not done : .FOF_Done=0.
*
*/
void fof_link_pseudo_heads(int target)
{
int prev,tail,head;
if (FOFDataGet[target].Task!=ThisTask)
{
printf("(%d) received ID=%d task=%d\n",ThisTask,FOFDataGet[target].ID,FOFDataGet[target].Task);
endrun(1956);
}
prev = FOFDataGet[target].FOF_Prev;
if (SphP[prev].FOF_Prev==-1) /* a dead end */
{
FOFDataResult[target].FOF_Prev = -1;
FOFDataResult[target].FOF_Head = -1;
FOFDataResult[target].FOF_Done = 1;
return;
}
/* find a valid head */
head = SphP[prev].FOF_Head;
if (head==-1)
endrun(1957);
if (SphP[head].FOF_Head==-1) /* a dead end */
{
FOFDataResult[target].FOF_Prev = -1;
FOFDataResult[target].FOF_Head = -1;
FOFDataResult[target].FOF_Done = 1;
return;
}
if (head!=SphP[head].FOF_Head)
{
printf("(%d) the head particle i=%d (i.FOF_Head=%d) is not a head !!! (i.FOF_CPUPrev=%d)\n",ThisTask,head,SphP[head].FOF_Head,SphP[head].FOF_CPUPrev);
endrun(1958);
}
/* head is now either
- a true head
- a pseudo head, but not done
- a pseudo head, done (prev=-2)
- a pseudo head, done
*/
/* a true head */
if ((SphP[head].FOF_CPUHead==ThisTask)&&(SphP[head].FOF_CPUPrev!=-2))
{
FOFDataResult[target].FOF_Prev = head; /* for this pseudo head, the head remains the same, but the real head is stored in prev */
FOFDataResult[target].FOF_CPUPrev = ThisTask;
FOFDataResult[target].FOF_Head = FOFDataGet[target].Index;
FOFDataResult[target].FOF_CPUHead = ThisTask;
FOFDataResult[target].FOF_Done = 1;
return;
}
/* a pseudo head, but not done */
if ( (SphP[head].FOF_CPUHead!=ThisTask) && (SphP[head].FOF_Done==0))
{
FOFDataResult[target].FOF_Prev = FOFDataGet[target].FOF_Prev; /* do nothing */
FOFDataResult[target].FOF_CPUPrev = ThisTask;
FOFDataResult[target].FOF_Head = FOFDataGet[target].FOF_Head; /* do nothing */
FOFDataResult[target].FOF_CPUHead = ThisTask;
FOFDataResult[target].FOF_Done = 0;
return;
}
/* a pseudo head, but done (prev=-2) */
if ((SphP[head].FOF_CPUHead==ThisTask)&&(SphP[head].FOF_CPUPrev==-2))
{
FOFDataResult[target].FOF_Prev = SphP[head].FOF_Prev;
FOFDataResult[target].FOF_CPUPrev = -2; /* this is to keep the trace that the particle was as pseudo head, as its new FOF_CPUHead could be now its own CPU */
FOFDataResult[target].FOF_Head = FOFDataGet[target].Index; /* WARNING : THOSE PARTICLES ARE PSEUDO HEAD THAT LINK TO OTHER CPU THAN THE INITIAL ONE */
FOFDataResult[target].FOF_CPUHead = SphP[head].FOF_CPUHead; /* instead of this ThisTask */
FOFDataResult[target].FOF_Done = 1;
return;
}
/* a pseudo head, but done */
if ( (SphP[head].FOF_CPUHead!=ThisTask) && (SphP[head].FOF_Done==1))
{
FOFDataResult[target].FOF_Prev = SphP[head].FOF_Prev;
FOFDataResult[target].FOF_CPUPrev = -2; /* this is to keep the trace that the particle was as pseudo head, as its new FOF_CPUHead could be now its own CPU */
FOFDataResult[target].FOF_Head = FOFDataGet[target].Index; /* WARNING : THOSE PARTICLES ARE PSEUDO HEAD THAT LINK TO OTHER CPU THAN THE INITIAL ONE */
FOFDataResult[target].FOF_CPUHead = SphP[head].FOF_CPUHead; /* instead of this ThisTask */
FOFDataResult[target].FOF_Done = 1;
return;
}
printf("we have a problem here !!\n");
endrun(1974);
}
/*! 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 fof_compare_key(const void *a, const void *b)
{
if(((struct FOFdata_in *) a)->Task < (((struct FOFdata_in *) b)->Task))
return -1;
if(((struct FOFdata_in *) a)->Task > (((struct FOFdata_in *) b)->Task))
return +1;
return 0;
}
/*! 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 fof_compare_phdata(const void *a, const void *b)
{
if(((struct fof_phdata *) a)->prev < (((struct fof_phdata *) b)->prev))
return -1;
if(((struct fof_phdata *) a)->prev > (((struct fof_phdata *) b)->prev))
return +1;
return 0;
}
/*! 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 fof_compare_groups_key(const void *a, const void *b)
{
if(((struct fof_groups_data *) a)->Task < (((struct fof_groups_data *) b)->Task))
return -1;
if(((struct fof_groups_data *) a)->Task > (((struct fof_groups_data *) b)->Task))
return +1;
return 0;
}
/*! This is a comparison kernel for a sort routine, which is used to sort
* particles according to some physical quantity.
*/
int fof_compare_particles(const void *a, const void *b)
{
/* if the mass is equal, sort with respect to ID */
if(((struct fof_particles_in *) a)->Mass == (((struct fof_particles_in *) b)->Mass))
{
if(((struct fof_particles_in *) a)->ID < (((struct fof_particles_in *) b)->ID))
return +1;
if(((struct fof_particles_in *) a)->ID > (((struct fof_particles_in *) b)->ID))
return -1;
return 0;
}
else
{
if(((struct fof_particles_in *) a)->Mass < (((struct fof_particles_in *) b)->Mass))
return +1;
if(((struct fof_particles_in *) a)->Mass > (((struct fof_particles_in *) b)->Mass))
return -1;
return 0;
}
}
/*! This function is used as a comparison kernel in a sort routine. It is
* used to group particles in the communication buffer that are going to
* be sent to the same CPU.
*/
int grav_fof_compare_key(const void *a, const void *b)
{
if(((struct FOFgravdata_index *) a)->Task < (((struct FOFgravdata_index *) b)->Task))
return -1;
if(((struct FOFgravdata_index *) a)->Task > (((struct FOFgravdata_index *) b)->Task))
return +1;
return 0;
}
#endif // FOF
c4science · Help