Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79502875
hydra.c
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Mon, Aug 26, 18:57
Size
20 KB
Mime Type
text/x-c
Expires
Wed, Aug 28, 18:57 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20197209
Attached To
rPNBODY pNbody
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
Log In to Comment