Homec4science
Diffusion Gear (master) Browse

Gear/src/stars_density.c
4b9ba9553668master

Last Change
Raw File

stars_density.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file density.c
* \brief SPH density computation and smoothing length determination
*
* This file contains the "first SPH loop", where the SPH densities and
* some auxiliary quantities are computed. If the number of neighbours
* obtained falls outside the target range, the correct smoothing
* length is determined iteratively, if needed.
*/
#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 CHIMIE
/*! This function computes the local density for each active SPH particle,
* the number of neighbours in the current smoothing radius, and the
* divergence and curl of the velocity field. The pressure is updated as
* well. If a particle with its smoothing region is fully inside the
* local domain, it is not exported to the other processors. The function
* also detects particles that have a number of neighbours outside the
* allowed tolerance range. For these particles, the smoothing length is
* adjusted accordingly, and the density computation is executed again.
* Note that the smoothing length is not allowed to fall below the lower
* bound set by MinGasHsml.
*/
void stars_density(void)
{
long long ntot, ntotleft;
int *noffset, *nbuffer, *nsend, *nsend_local, *numlist, *ndonelist;
int i, j, n,m, ndone, npleft, maxfill, source, iter = 0;
int level, ngrp, sendTask, recvTask, place, nexport;
double tstart, tend, tstart_ngb = 0, tend_ngb = 0;
double sumt, sumcomm, timengb, sumtimengb;
double timecomp = 0, timeimbalance = 0, timecommsumm = 0, sumimbalance;
MPI_Status status;
#ifdef DETAILED_CPU_OUTPUT_IN_STARS_DENSITY
double *timengblist;
double *timecomplist;
double *timecommsummlist;
double *timeimbalancelist;
#endif
#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
if (ThisTask==0)
printf("Start star density computation.\n");
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);
/* `NumStUpdate' gives the number of particles on this processor that want a chimie computation */
for(n = N_gas, NumStUpdate = 0; n < N_gas+N_stars; n++)
{
m = P[n].StPIdx;
StP[m].Left = StP[m].Right = 0;
#ifdef AVOIDNUMNGBPROBLEM
StP[m].OldNumNgb = -1;
#endif
if(P[n].Ti_endstep == All.Ti_Current)
NumStUpdate++;
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
/* we will repeat the whole thing for those particles where we didn't
* find enough neighbours
*/
do
{
i = N_gas; /* beginn with the first star */
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 */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas+N_stars && nexport < All.BunchSizeDensity - NTask; i++)
if(P[i].Ti_endstep == All.Ti_Current)
{
if(P[i].Type != ST)
{
printf("P[i].Type != ST, we better stop.\n");
printf("N_gas=%d (type=%d) i=%d (type=%d)\n",N_gas,P[N_gas].Type,i,P[i].Type);
printf("Please, check that you do not use PEANOHILBERT\n");
endrun(666001);
}
ndone++;
m = P[i].StPIdx;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
stars_density_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
StarsDensDataIn[nexport].Pos[0] = P[i].Pos[0];
StarsDensDataIn[nexport].Pos[1] = P[i].Pos[1];
StarsDensDataIn[nexport].Pos[2] = P[i].Pos[2];
StarsDensDataIn[nexport].Hsml = StP[m].Hsml;
StarsDensDataIn[nexport].Index = i;
StarsDensDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
}
}
tend = second();
timecomp += timediff(tstart, tend);
qsort(StarsDensDataIn, nexport, sizeof(struct starsdensdata_in), stars_dens_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.BunchSizeDensity)
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(&StarsDensDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct starsdensdata_in), MPI_BYTE,
recvTask, TAG_DENS_A,
&StarsDensDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct starsdensdata_in),
MPI_BYTE, recvTask, TAG_DENS_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);
tstart = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
stars_density_evaluate(j, 1);
tend = second();
timecomp += timediff(tstart, tend);
/* do a block to explicitly 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.BunchSizeDensity)
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(&StarsDensDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct starsdensdata_out),
MPI_BYTE, recvTask, TAG_DENS_B,
&StarsDensDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct starsdensdata_out),
MPI_BYTE, recvTask, TAG_DENS_B, MPI_COMM_WORLD, &status);
/* add the result to the particles */
for(j = 0; j < nsend_local[recvTask]; j++)
{
source = j + noffset[recvTask];
place = P[StarsDensDataIn[source].Index].StPIdx;
StP[place].NumNgb += StarsDensDataPartialResult[source].Ngb;
StP[place].Density += StarsDensDataPartialResult[source].Rho;
StP[place].Volume += StarsDensDataPartialResult[source].Volume;
#ifdef CHIMIE_KINETIC_FEEDBACK
StP[place].NgbMass += StarsDensDataPartialResult[source].NgbMass;
#endif
StP[place].DhsmlDensityFactor += StarsDensDataPartialResult[source].DhsmlDensity;
}
}
}
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];
}
/* do final operations on results */
tstart = second();
for(i = N_gas, npleft = 0; i < N_gas+N_stars; i++)
{
if((P[i].Ti_endstep == All.Ti_Current))
{
m = P[i].StPIdx;
StP[m].DhsmlDensityFactor =
1 / (1 + StP[m].Hsml * StP[m].DhsmlDensityFactor / (NUMDIMS * StP[m].Density));
/* now check whether we had enough neighbours */
if(StP[m].NumNgb < (All.DesNumNgb - All.MaxNumNgbDeviation) ||
(StP[m].NumNgb > (All.DesNumNgb + All.MaxNumNgbDeviation)
&& StP[m].Hsml > (1.01 * All.MinGasHsml)))
{
#ifdef AVOIDNUMNGBPROBLEM
// if((StP[m].NumNgb>StP[m].OldNumNgb-All.MaxNumNgbDeviation/10000.)
// &&(StP[m].NumNgb<StP[m].OldNumNgb+All.MaxNumNgbDeviation/10000.))
if(StP[m].NumNgb==StP[m].OldNumNgb)
{
P[i].Ti_endstep = -P[i].Ti_endstep - 1;
printf("ID=%d NumNgb=%g OldNumNgb=%g\n",P[i].ID,StP[m].NumNgb,StP[m].OldNumNgb);
continue;
}
StP[m].OldNumNgb = StP[m].NumNgb;
#endif
/* need to redo this particle */
npleft++;
if(StP[m].Left > 0 && StP[m].Right > 0)
if((StP[m].Right - StP[m].Left) < 1.0e-3 * StP[m].Left)
{
/* this one should be ok */
npleft--;
P[i].Ti_endstep = -P[i].Ti_endstep - 1; /* Mark as inactive */
continue;
}
if(StP[m].NumNgb < (All.DesNumNgb - All.MaxNumNgbDeviation))
StP[m].Left = dmax(StP[m].Hsml, StP[m].Left);
else
{
if(StP[m].Right != 0)
{
if(StP[m].Hsml < StP[m].Right)
StP[m].Right = StP[m].Hsml;
}
else
StP[m].Right = StP[m].Hsml;
}
if(iter >= MAXITER - 10)
{
printf
("i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n pos=(%g|%g|%g)\n",
i, ThisTask, (int) P[i].ID, StP[m].Hsml, StP[m].Left, StP[m].Right,
(float) StP[m].NumNgb, StP[m].Right - StP[m].Left, P[i].Pos[0], P[i].Pos[1],
P[i].Pos[2]);
fflush(stdout);
}
if(StP[m].Right > 0 && StP[m].Left > 0)
StP[m].Hsml = pow(0.5 * (pow(StP[m].Left, 3) + pow(StP[m].Right, 3)), 1.0 / 3);
else
{
if(StP[m].Right == 0 && StP[m].Left == 0)
{
printf
("i=%d task=%d ID=%d Hsml=%g Left=%g Right=%g Ngbs=%g Right-Left=%g\n pos=(%g|%g|%g)\n",
i, ThisTask, (int) P[i].ID, StP[m].Hsml, StP[m].Left, StP[m].Right,
(float) StP[m].NumNgb, StP[m].Right - StP[m].Left, P[i].Pos[0], P[i].Pos[1],
P[i].Pos[2]);
fflush(stdout);
endrun(8188); /* can't occur */
}
if(StP[m].Right == 0 && StP[m].Left > 0)
{
if(P[i].Type == ST && fabs(StP[m].NumNgb - All.DesNumNgb) < 0.5 * All.DesNumNgb)
{
StP[m].Hsml *=
1 - (StP[m].NumNgb -
All.DesNumNgb) / (NUMDIMS * StP[m].NumNgb) * StP[m].DhsmlDensityFactor;
}
else
StP[m].Hsml *= 1.26;
}
if(StP[m].Right > 0 && StP[m].Left == 0)
{
if(P[i].Type == ST && fabs(StP[m].NumNgb - All.DesNumNgb) < 0.5 * All.DesNumNgb)
{
StP[m].Hsml *=
1 - (StP[m].NumNgb -
All.DesNumNgb) / (NUMDIMS * StP[m].NumNgb) * StP[m].DhsmlDensityFactor;
}
else
StP[m].Hsml /= 1.26;
}
}
if(StP[m].Hsml < All.MinGasHsml)
StP[m].Hsml = All.MinGasHsml;
}
else
P[i].Ti_endstep = -P[i].Ti_endstep - 1; /* Mark as inactive */
}
}
tend = second();
timecomp += timediff(tstart, tend);
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)
{
if(iter == 0)
tstart_ngb = second();
iter++;
if(iter > 0 && ThisTask == 0)
{
printf("ngb iteration %d: need to repeat for %d%09d particles.\n", iter,
(int) (ntot / 1000000000), (int) (ntot % 1000000000));
fflush(stdout);
}
if(iter > MAXITER)
{
printf("failed to converge in neighbour iteration in density()\n");
fflush(stdout);
endrun(1155);
}
}
else
tend_ngb = second();
}
while(ntot > 0);
/* mark as active again */
for(i = 0; i < NumPart; i++)
if(P[i].Ti_endstep < 0)
P[i].Ti_endstep = -P[i].Ti_endstep - 1;
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
/* collect some timing information */
if(iter > 0)
timengb = timediff(tstart_ngb, tend_ngb);
else
timengb = 0;
MPI_Reduce(&timengb, &sumtimengb, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timecomp, &sumt, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timecommsumm, &sumcomm, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&timeimbalance, &sumimbalance, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#ifdef DETAILED_CPU_OUTPUT_IN_STARS_DENSITY
numlist = malloc(sizeof(int) * NTask);
timengblist = malloc(sizeof(double) * NTask);
timecomplist = malloc(sizeof(double) * NTask);
timecommsummlist = malloc(sizeof(double) * NTask);
timeimbalancelist = malloc(sizeof(double) * NTask);
MPI_Gather(&NumStUpdate, 1, MPI_INT, numlist, 1, MPI_INT, 0, MPI_COMM_WORLD);
MPI_Gather(&timengb, 1, MPI_DOUBLE, timengblist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(&timecomp, 1, MPI_DOUBLE, timecomplist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(&timecommsumm, 1, MPI_DOUBLE, timecommsummlist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
MPI_Gather(&timeimbalance, 1, MPI_DOUBLE, timeimbalancelist, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
if(ThisTask == 0)
{
fprintf(FdTimings, "\n stars density \n\n");
fprintf(FdTimings, "Nupdate ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12d ",numlist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timengb ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timengblist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timecomp ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timecomplist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timecommsumm ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timecommsummlist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "timeimbalance ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12g ",timeimbalancelist[i]);
fprintf(FdTimings, "\n");
fprintf(FdTimings, "\n");
}
free(timeimbalancelist);
free(timecommsummlist);
free(timecomplist);
free(timengblist);
free(numlist);
#endif
if(ThisTask == 0)
{
All.CPU_ChimieDensCompWalk += sumt / NTask;
All.CPU_ChimieDensCommSumm += sumcomm / NTask;
All.CPU_ChimieDensImbalance += sumimbalance / NTask;
All.CPU_ChimieDensEnsureNgb += sumtimengb / NTask;
}
if (ThisTask==0)
printf("Stars density computation done.\n");
}
/*! This function represents the core of the SPH density computation. The
* target particle may either be local, or reside in the communication
* buffer.
*/
void stars_density_evaluate(int target, int mode)
{
int j, n, startnode, numngb, numngb_inbox;
double h, h2, hinv, hinv3, hinv4;
double rho, volume, wk, dwk;
double dx, dy, dz, r, r2, u, mass_j,rho_j;
double weighted_numngb, dhsmlrho;
int target_stp;
FLOAT *pos;
#ifdef CHIMIE_KINETIC_FEEDBACK
double ngbmass;
#endif
if(mode == 0)
{
pos = P[target].Pos;
target_stp = P[target].StPIdx;
h = StP[target_stp].Hsml;
}
else
{
pos = StarsDensDataGet[target].Pos;
h = StarsDensDataGet[target].Hsml;
}
h2 = h * h;
hinv = 1.0 / h;
#ifndef TWODIMS
hinv3 = hinv * hinv * hinv;
#else
hinv3 = hinv * hinv / boxSize_Z;
#endif
hinv4 = hinv3 * hinv;
rho = 0;
volume = 0;
weighted_numngb = 0;
dhsmlrho = 0;
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass=0;
#endif
startnode = All.MaxPart;
numngb = 0;
do
{
numngb_inbox = ngb_treefind_variable_for_chimie(&pos[0], h, &startnode);
for(n = 0; n < numngb_inbox; 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 /* now 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)
{
numngb++;
r = sqrt(r2);
u = r * hinv;
if(u < 0.5)
{
wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
dwk = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
}
else
{
wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
dwk = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
}
mass_j = P[j].Mass;
rho_j = SphP[j].Density;
rho += mass_j * wk;
volume += mass_j/rho_j * wk;
weighted_numngb += NORM_COEFF * wk / hinv3;
dhsmlrho += -mass_j * (NUMDIMS * hinv * wk + u * dwk);
#ifdef CHIMIE_KINETIC_FEEDBACK
ngbmass += mass_j;
#endif
}
}
}
while(startnode >= 0);
if(mode == 0)
{
StP[target_stp].NumNgb = weighted_numngb;
StP[target_stp].Density = rho;
StP[target_stp].Volume = volume;
StP[target_stp].DhsmlDensityFactor = dhsmlrho;
#ifdef CHIMIE_KINETIC_FEEDBACK
StP[target_stp].NgbMass = ngbmass;
#endif
}
else
{
StarsDensDataResult[target].Rho = rho;
StarsDensDataResult[target].Volume = volume;
StarsDensDataResult[target].Ngb = weighted_numngb;
StarsDensDataResult[target].DhsmlDensity = dhsmlrho;
#ifdef CHIMIE_KINETIC_FEEDBACK
StarsDensDataResult[target].NgbMass = ngbmass;
#endif
}
}
/*! This routine is a comparison kernel used in a sort routine to group
* particles that are exported to the same processor.
*/
int stars_dens_compare_key(const void *a, const void *b)
{
if(((struct starsdensdata_in *) a)->Task < (((struct starsdensdata_in *) b)->Task))
return -1;
if(((struct starsdensdata_in *) a)->Task > (((struct starsdensdata_in *) b)->Task))
return +1;
return 0;
}
#endif
c4science · Help