Page MenuHomec4science
Files F79502875

hydra.c
No OneTemporary
Actions

File Metadata

Created
Mon, Aug 26, 18:57

hydra.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"
/*! \file hydra.c
* \brief Computation of SPH forces and rate of entropy generation
*
* This file contains the "second SPH loop", where the SPH forces are
* computed, and where the rate of change of entropy due to the shock heating
* (via artificial viscosity) is computed.
*/
static double hubble_a, atime, hubble_a2, fac_mu, fac_vsic_fix, a3inv, fac_egy;
#ifdef FEEDBACK
static double fac_pow;
#endif
#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
/*! This function is the driver routine for the calculation of hydrodynamical
* force and rate of change of entropy due to shock heating for all active
* particles .
*/
void hydro_force(void)
{
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 soundspeed_i;
double tstart, tend, sumt, sumcomm;
double timecomp = 0, timecommsumm = 0, timeimbalance = 0, sumimbalance;
MPI_Status status;
#ifdef PERIODIC
boxSize = All.BoxSize;
boxHalf = 0.5 * All.BoxSize;
#ifdef LONG_X
boxHalf_X = boxHalf * LONG_X;
boxSize_X = boxSize * LONG_X;
#endif
#ifdef LONG_Y
boxHalf_Y = boxHalf * LONG_Y;
boxSize_Y = boxSize * LONG_Y;
#endif
#ifdef LONG_Z
boxHalf_Z = boxHalf * LONG_Z;
boxSize_Z = boxSize * LONG_Z;
#endif
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
double v1m,v2m;
#endif
if(All.ComovingIntegrationOn)
{
/* Factors for comoving integration of hydro */
hubble_a = All.Omega0 / (All.Time * All.Time * All.Time)
+ (1 - All.Omega0 - All.OmegaLambda) / (All.Time * All.Time) + All.OmegaLambda;
hubble_a = All.Hubble * sqrt(hubble_a);
hubble_a2 = All.Time * All.Time * hubble_a;
fac_mu = pow(All.Time, 3 * (GAMMA - 1) / 2) / All.Time;
fac_egy = pow(All.Time, 3 * (GAMMA - 1));
fac_vsic_fix = hubble_a * pow(All.Time, 3 * GAMMA_MINUS1);
a3inv = 1 / (All.Time * All.Time * All.Time);
atime = All.Time;
#ifdef FEEDBACK
fac_pow = fac_egy*atime*atime;
#endif
}
else
{
hubble_a = hubble_a2 = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0;
#ifdef FEEDBACK
fac_pow = 1.0;
#endif
}
/* `NumSphUpdate' gives the number of particles on this processor that want a force update */
for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
{
#ifdef SFR
if((P[n].Ti_endstep == All.Ti_Current) && (P[n].Type == 0))
#else
if(P[n].Ti_endstep == All.Ti_Current)
#endif
#ifdef MULTIPHASE
if(SphP[n].Phase == GAS_SPH)
#endif
NumSphUpdate++;
}
numlist = malloc(NTask * sizeof(int) * NTask);
MPI_Allgather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, MPI_COMM_WORLD);
for(i = 0, ntot = 0; i < NTask; i++)
ntot += numlist[i];
free(numlist);
noffset = malloc(sizeof(int) * NTask); /* offsets of bunches in common list */
nbuffer = malloc(sizeof(int) * NTask);
nsend_local = malloc(sizeof(int) * NTask);
nsend = malloc(sizeof(int) * NTask * NTask);
ndonelist = malloc(sizeof(int) * NTask);
i = 0; /* first particle for this task */
ntotleft = ntot; /* particles left for all tasks together */
while(ntotleft > 0)
{
for(j = 0; j < NTask; j++)
nsend_local[j] = 0;
/* do local particles and prepare export list */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeHydro - NTask; i++)
#ifdef SFR
if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
if(P[i].Ti_endstep == All.Ti_Current)
#endif
{
#ifdef MULTIPHASE
if(SphP[i].Phase == GAS_SPH)
{
#endif
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
hydro_evaluate(i, 0);
for(j = 0; j < NTask; j++)
{
if(Exportflag[j])
{
for(k = 0; k < 3; k++)
{
HydroDataIn[nexport].Pos[k] = P[i].Pos[k];
HydroDataIn[nexport].Vel[k] = SphP[i].VelPred[k];
}
HydroDataIn[nexport].Hsml = SphP[i].Hsml;
#ifdef FEEDBACK
HydroDataIn[nexport].EnergySN = SphP[i].EnergySN;
#endif
HydroDataIn[nexport].Mass = P[i].Mass;
HydroDataIn[nexport].DhsmlDensityFactor = SphP[i].DhsmlDensityFactor;
HydroDataIn[nexport].Density = SphP[i].Density;
HydroDataIn[nexport].Pressure = SphP[i].Pressure;
HydroDataIn[nexport].Timestep = P[i].Ti_endstep - P[i].Ti_begstep;
#ifdef WITH_ID_IN_HYDRA
HydroDataIn[nexport].ID = P[i].ID;
#endif
/* calculation of F1 */
soundspeed_i = sqrt(GAMMA * SphP[i].Pressure / SphP[i].Density);
HydroDataIn[nexport].F1 = fabs(SphP[i].DivVel) /
(fabs(SphP[i].DivVel) + SphP[i].CurlVel +
0.0001 * soundspeed_i / SphP[i].Hsml / fac_mu);
HydroDataIn[nexport].Index = i;
HydroDataIn[nexport].Task = j;
nexport++;
nsend_local[j]++;
}
}
#ifdef MULTIPHASE
}
#endif
}
tend = second();
timecomp += timediff(tstart, tend);
qsort(HydroDataIn, nexport, sizeof(struct hydrodata_in), hydro_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.BunchSizeHydro)
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(&HydroDataIn[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct hydrodata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A,
&HydroDataGet[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct hydrodata_in), MPI_BYTE,
recvTask, TAG_HYDRO_A, MPI_COMM_WORLD, &status);
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
/* now do the imported particles */
tstart = second();
for(j = 0; j < nbuffer[ThisTask]; j++)
hydro_evaluate(j, 1);
tend = second();
timecomp += timediff(tstart, tend);
/* do a block to measure imbalance */
tstart = second();
MPI_Barrier(MPI_COMM_WORLD);
tend = second();
timeimbalance += timediff(tstart, tend);
/* get the result */
tstart = second();
for(j = 0; j < NTask; j++)
nbuffer[j] = 0;
for(ngrp = level; ngrp < (1 << PTask); ngrp++)
{
maxfill = 0;
for(j = 0; j < NTask; j++)
{
if((j ^ ngrp) < NTask)
if(maxfill < nbuffer[j] + nsend[(j ^ ngrp) * NTask + j])
maxfill = nbuffer[j] + nsend[(j ^ ngrp) * NTask + j];
}
if(maxfill >= All.BunchSizeHydro)
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(&HydroDataResult[nbuffer[ThisTask]],
nsend[recvTask * NTask + ThisTask] * sizeof(struct hydrodata_out),
MPI_BYTE, recvTask, TAG_HYDRO_B,
&HydroDataPartialResult[noffset[recvTask]],
nsend_local[recvTask] * sizeof(struct hydrodata_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 = HydroDataIn[source].Index;
for(k = 0; k < 3; k++)
SphP[place].HydroAccel[k] += HydroDataPartialResult[source].Acc[k];
SphP[place].DtEntropy += HydroDataPartialResult[source].DtEntropy;
#ifdef FEEDBACK
SphP[place].DtEgySpecFeedback += HydroDataPartialResult[source].DtEgySpecFeedback;
#endif
if(SphP[place].MaxSignalVel < HydroDataPartialResult[source].MaxSignalVel)
SphP[place].MaxSignalVel = HydroDataPartialResult[source].MaxSignalVel;
#ifdef COMPUTE_VELOCITY_DISPERSION
for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
SphP[place].VelocityDispersion[k] += HydroDataPartialResult[source].VelocityDispersion[k];
#endif
}
}
}
for(j = 0; j < NTask; j++)
if((j ^ ngrp) < NTask)
nbuffer[j] += nsend[(j ^ ngrp) * NTask + j];
}
tend = second();
timecommsumm += timediff(tstart, tend);
level = ngrp - 1;
}
MPI_Allgather(&ndone, 1, MPI_INT, ndonelist, 1, MPI_INT, MPI_COMM_WORLD);
for(j = 0; j < NTask; j++)
ntotleft -= ndonelist[j];
}
free(ndonelist);
free(nsend);
free(nsend_local);
free(nbuffer);
free(noffset);
/* do final operations on results */
tstart = second();
for(i = 0; i < N_gas; i++)
#ifdef SFR
if((P[i].Ti_endstep == All.Ti_Current) && (P[i].Type == 0))
#else
if(P[i].Ti_endstep == All.Ti_Current)
#endif
{
SphP[i].DtEntropy *= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].Density, GAMMA_MINUS1));
#ifdef SPH_BND_PARTICLES
if(P[i].ID == 0)
{
SphP[i].DtEntropy = 0;
for(k = 0; k < 3; k++)
SphP[i].HydroAccel[k] = 0;
}
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
if (SphP[i].VelocityDispersion[0] != 0)
{
/* compute sigma r */
v1m = SphP[i].VelocityDispersion[1]/SphP[i].VelocityDispersion[0];
v2m = SphP[i].VelocityDispersion[2]/SphP[i].VelocityDispersion[0];
if (v2m > v1m*v1m)
SphP[i].OptVar = sqrt(v2m - v1m*v1m);
else
SphP[i].OptVar = 0.0;
}
else
SphP[i].OptVar = 0.0;
#endif
}
tend = second();
timecomp += timediff(tstart, tend);
/* collect some timing information */
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);
if(ThisTask == 0)
{
All.CPU_HydCompWalk += sumt / NTask;
All.CPU_HydCommSumm += sumcomm / NTask;
All.CPU_HydImbalance += sumimbalance / NTask;
}
}
/*! This function is the 'core' of the SPH force computation. A target
* particle is specified which may either be local, or reside in the
* communication buffer.
*/
void hydro_evaluate(int target, int mode)
{
int j, k, n, timestep, startnode, numngb;
FLOAT *pos, *vel;
FLOAT mass, h_i, dhsmlDensityFactor, rho, pressure, f1, f2;
double acc[3], dtEntropy, maxSignalVel;
double dx, dy, dz, dvx, dvy, dvz;
double h_i2, hinv=1, hinv4;
double p_over_rho2_i, p_over_rho2_j, soundspeed_i, soundspeed_j;
double hfc, dwk_i, vdotr, vdotr2, visc, mu_ij, rho_ij=0, vsig;
double h_j, dwk_j, r, r2, u=0, hfc_visc;
int phase=0;
#ifdef FEEDBACK
int EnergySN;
double wk,wk_i,wk_j,uintspec,hinv3;
double dtEgySpecFeedback=0;
#endif
#ifdef COMPUTE_VELOCITY_DISPERSION
double VelocityDispersion[VELOCITY_DISPERSION_SIZE];
for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
VelocityDispersion[k]=0.0;
#endif
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
if(mode == 0)
{
pos = P[target].Pos;
vel = SphP[target].VelPred;
h_i = SphP[target].Hsml;
#ifdef FEEDBACK
EnergySN = SphP[target].EnergySN;
#endif
mass = P[target].Mass;
dhsmlDensityFactor = SphP[target].DhsmlDensityFactor;
rho = SphP[target].Density;
pressure = SphP[target].Pressure;
timestep = P[target].Ti_endstep - P[target].Ti_begstep;
soundspeed_i = sqrt(GAMMA * pressure / rho);
f1 = fabs(SphP[target].DivVel) /
(fabs(SphP[target].DivVel) + SphP[target].CurlVel +
0.0001 * soundspeed_i / SphP[target].Hsml / fac_mu);
}
else
{
pos = HydroDataGet[target].Pos;
vel = HydroDataGet[target].Vel;
h_i = HydroDataGet[target].Hsml;
#ifdef FEEDBACK
EnergySN = HydroDataGet[target].EnergySN;
#endif
mass = HydroDataGet[target].Mass;
dhsmlDensityFactor = HydroDataGet[target].DhsmlDensityFactor;
rho = HydroDataGet[target].Density;
pressure = HydroDataGet[target].Pressure;
timestep = HydroDataGet[target].Timestep;
soundspeed_i = sqrt(GAMMA * pressure / rho);
f1 = HydroDataGet[target].F1;
}
/* initialize variables before SPH loop is started */
acc[0] = acc[1] = acc[2] = dtEntropy = 0;
maxSignalVel = 0;
#ifdef FEEDBACK
dtEgySpecFeedback=0;
#endif
p_over_rho2_i = pressure / (rho * rho) * dhsmlDensityFactor;
h_i2 = h_i * h_i;
/* Now start the actual SPH computation for this particle */
startnode = All.MaxPart;
do
{
numngb = ngb_treefind_pairs(&pos[0], h_i, phase, &startnode);
for(n = 0; n < numngb; n++)
{
j = Ngblist[n];
dx = pos[0] - P[j].Pos[0];
dy = pos[1] - P[j].Pos[1];
dz = pos[2] - P[j].Pos[2];
#ifdef PERIODIC /* find the closest image in the given box size */
if(dx > boxHalf_X)
dx -= boxSize_X;
if(dx < -boxHalf_X)
dx += boxSize_X;
if(dy > boxHalf_Y)
dy -= boxSize_Y;
if(dy < -boxHalf_Y)
dy += boxSize_Y;
if(dz > boxHalf_Z)
dz -= boxSize_Z;
if(dz < -boxHalf_Z)
dz += boxSize_Z;
#endif
r2 = dx * dx + dy * dy + dz * dz;
h_j = SphP[j].Hsml;
if(r2 < h_i2 || r2 < h_j * h_j)
{
r = sqrt(r2);
if(r > 0)
{
p_over_rho2_j = SphP[j].Pressure / (SphP[j].Density * SphP[j].Density);
soundspeed_j = sqrt(GAMMA * p_over_rho2_j * SphP[j].Density);
dvx = vel[0] - SphP[j].VelPred[0];
dvy = vel[1] - SphP[j].VelPred[1];
dvz = vel[2] - SphP[j].VelPred[2];
vdotr = dx * dvx + dy * dvy + dz * dvz;
if(All.ComovingIntegrationOn)
vdotr2 = vdotr + hubble_a2 * r2;
else
vdotr2 = vdotr;
if(r2 < h_i2)
{
hinv = 1.0 / h_i;
#ifndef TWODIMS
hinv4 = hinv * hinv * hinv * hinv;
#else
hinv4 = hinv * hinv * hinv / boxSize_Z;
#endif
u = r * hinv;
if(u < 0.5)
dwk_i = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
else
dwk_i = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
}
else
{
dwk_i = 0;
}
if(r2 < h_j * h_j)
{
hinv = 1.0 / h_j;
#ifndef TWODIMS
hinv4 = hinv * hinv * hinv * hinv;
#else
hinv4 = hinv * hinv * hinv / boxSize_Z;
#endif
u = r * hinv;
if(u < 0.5)
dwk_j = hinv4 * u * (KERNEL_COEFF_3 * u - KERNEL_COEFF_4);
else
dwk_j = hinv4 * KERNEL_COEFF_6 * (1.0 - u) * (1.0 - u);
}
else
{
dwk_j = 0;
}
if(soundspeed_i + soundspeed_j > maxSignalVel)
maxSignalVel = soundspeed_i + soundspeed_j;
if(vdotr2 < 0) /* ... artificial viscosity */
{
mu_ij = fac_mu * vdotr2 / r; /* note: this is negative! */
vsig = soundspeed_i + soundspeed_j - 3 * mu_ij;
if(vsig > maxSignalVel)
maxSignalVel = vsig;
rho_ij = 0.5 * (rho + SphP[j].Density);
f2 =
fabs(SphP[j].DivVel) / (fabs(SphP[j].DivVel) + SphP[j].CurlVel +
0.0001 * soundspeed_j / fac_mu / SphP[j].Hsml);
visc = 0.25 * All.ArtBulkViscConst * vsig * (-mu_ij) / rho_ij * (f1 + f2);
/* .... end artificial viscosity evaluation */
#ifndef NOVISCOSITYLIMITER
/* make sure that viscous acceleration is not too large */
dt = imax(timestep, (P[j].Ti_endstep - P[j].Ti_begstep)) * All.Timebase_interval;
if(dt > 0 && (dwk_i + dwk_j) < 0)
{
visc = dmin(visc, 0.5 * fac_vsic_fix * vdotr2 /
(0.5 * (mass + P[j].Mass) * (dwk_i + dwk_j) * r * dt));
}
#endif
}
else
visc = 0;
p_over_rho2_j *= SphP[j].DhsmlDensityFactor;
hfc_visc = 0.5 * P[j].Mass * visc * (dwk_i + dwk_j) / r;
hfc = hfc_visc + P[j].Mass * (p_over_rho2_i * dwk_i + p_over_rho2_j * dwk_j) / r;
acc[0] -= hfc * dx;
acc[1] -= hfc * dy;
acc[2] -= hfc * dz;
dtEntropy += 0.5 * hfc_visc * vdotr2;
/*****************************************/
/* FEEDBACK INTERACTION */
/*****************************************/
#ifdef FEEDBACK
rho_ij = 0.5 * (rho + SphP[j].Density);
if(P[j].Ti_endstep == All.Ti_Current)
{
/* additional feedback entropy */
if ((EnergySN > 0)||(SphP[j].EnergySN > 0))
{
/* find the thermal specific energy to release */
uintspec = 0.;
if (EnergySN > 0)
{
uintspec = 0;
}
if (SphP[j].EnergySN > 0)
{
uintspec += SphP[j].EnergySN * (1-All.SupernovaFractionInEgyKin)/ mass;
}
if(r2 < h_i2)
{
hinv = 1.0 / h_i;
hinv3 = hinv * hinv * hinv ;
u = r * hinv;
if(u < 0.5)
wk_i = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
else
wk_i = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
}
else
wk_i = 0;
if(r2 < h_j * h_j)
{
hinv = 1.0 / h_j;
hinv3 = hinv * hinv * hinv ;
u = r * hinv;
if(u < 0.5)
wk_j = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u);
else
wk_j = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u);
}
else
wk_j = 0;
wk = 0.5*(wk_i+wk_j);
/* dt in physical units */
dt = imax(timestep, (P[j].Ti_endstep - P[j].Ti_begstep)) * All.Timebase_interval/ hubble_a;
if(dt <= 0);
uintspec = 0;
/* thermal feedback */
uintspec = (uintspec/dt) *wk *0.5*(mass + P[j].Mass) /rho_ij *fac_pow ;
dtEntropy += uintspec;
dtEgySpecFeedback += uintspec;
}
}
#endif /* FEEDBACK */
}
}
}
}
while(startnode >= 0);
/* Now collect the result at the right place */
if(mode == 0)
{
for(k = 0; k < 3; k++)
SphP[target].HydroAccel[k] = acc[k];
SphP[target].DtEntropy = dtEntropy;
#ifdef FEEDBACK
SphP[target].DtEgySpecFeedback = dtEgySpecFeedback;
#endif
SphP[target].MaxSignalVel = maxSignalVel;
#ifdef COMPUTE_VELOCITY_DISPERSION
for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
SphP[target].VelocityDispersion[k] = VelocityDispersion[k];
#endif
}
else
{
for(k = 0; k < 3; k++)
HydroDataResult[target].Acc[k] = acc[k];
HydroDataResult[target].DtEntropy = dtEntropy;
#ifdef FEEDBACK
HydroDataResult[target].DtEgySpecFeedback = dtEgySpecFeedback;
#endif
HydroDataResult[target].MaxSignalVel = maxSignalVel;
#ifdef COMPUTE_VELOCITY_DISPERSION
for(k = 0; k < VELOCITY_DISPERSION_SIZE; k++)
HydroDataResult[target].VelocityDispersion[k] = VelocityDispersion[k];
#endif
}
}
/*! 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 hydro_compare_key(const void *a, const void *b)
{
if(((struct hydrodata_in *) a)->Task < (((struct hydrodata_in *) b)->Task))
return -1;
if(((struct hydrodata_in *) a)->Task > (((struct hydrodata_in *) b)->Task))
return +1;
return 0;
}

Event Timeline

c4science · Help