stars_density.c
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Aug 17, 11:52

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

stars_density.cNo OneTemporaryActions

File Metadata

stars_density.cView Options

Event Timeline

stars_density.c
No OneTemporary
Actions

stars_density.c
View Options