Homec4science
Diffusion Gear (master) Browse

Gear/src/sigvel.c
5ec004b6ec8dmaster

Last Change
Raw File

sigvel.c
View Options

/* ################################################################################## */
/* ### ### */
/* ### sigvel.c ### */
/* ### ### */
/* ### Original: hydra.c (public version of Gadget 2) ### */
/* ### Author: Volker Springel ### */
/* ### ### */
/* ### Modified: February 2011 ### */
/* ### Author: Fabrice Durier ### */
/* ### ### */
/* ################################################################################## */
#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 sigvel.c
* Re-Computation of SPH acceleration and maximum signal velocity
* for 'feedback' particles only.
*
* This file is a modified version of hydra.c, where the SPH forces are
* computed, and where the rate of change of entropy due to the shock heating
* (via artificial viscosity) is computed.
*/
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
#ifdef CHIMIE_THERMAL_FEEDBACK
static double hubble_a, atime, hubble_a2, fac_mu, fac_vsic_fix, a3inv, fac_egy;
#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 update of the calculation
* of hydrodynamical force due to the feedback impact
* from either SNe or BHs on active gas particles.
*/
void get_sigvel(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
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;
}
else
hubble_a = hubble_a2 = atime = fac_mu = fac_vsic_fix = a3inv = fac_egy = 1.0;
/* `NumSphUpdate' gives the number of feedback particles on this processor that want a force update */
for(n = 0, NumSphUpdate = 0; n < N_gas; n++)
{
if(P[n].Ti_endstep == All.Ti_Current && SphP[n].DeltaEgySpec > 0)
NumSphUpdate++;
for(j = 0; j < 3; j++)
SphP[n].FeedbackUpdatedAccel[j] = 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 && ntot)
printf("\t ---> Number of feedback particles = %ld \n\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 */
tstart = second();
for(nexport = 0, ndone = 0; i < N_gas && nexport < All.BunchSizeHydro - NTask; i++)
if(P[i].Ti_endstep == All.Ti_Current && SphP[i].DeltaEgySpec > 0)
{
ndone++;
for(j = 0; j < NTask; j++)
Exportflag[j] = 0;
get_sigvel_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;
HydroDataIn[nexport].Mass = P[i].Mass;
#ifdef DENSITY_INDEPENDENT_SPH
HydroDataIn[nexport].EgyRho = SphP[i].EgyWtDensity;
HydroDataIn[nexport].EntVarPred = SphP[i].EntVarPred;
HydroDataIn[nexport].DhsmlDensityFactor = SphP[i].DhsmlEgyDensityFactor;
#else
HydroDataIn[nexport].DhsmlDensityFactor = SphP[i].DhsmlDensityFactor;
#endif
HydroDataIn[nexport].Density = SphP[i].Density;
#ifdef DENSITY_INDEPENDENT_SPH
HydroDataIn[nexport].Pressure = updated_pressure(SphP[i].EntropyPred,SphP[i].EgyWtDensity,SphP[i].DeltaEgySpec);
#else
HydroDataIn[nexport].Pressure = updated_pressure(SphP[i].EntropyPred,SphP[i].Density,SphP[i].DeltaEgySpec);
#endif
HydroDataIn[nexport].Timestep = P[i].Ti_endstep - P[i].Ti_begstep;
/* calculation of F1 */
soundspeed_i = sqrt(GAMMA * updated_pressure(SphP[i].EntropyPred,SphP[i].Density,SphP[i].DeltaEgySpec) / 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]++;
}
}
}
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++)
get_sigvel_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].FeedbackUpdatedAccel[k] += HydroDataPartialResult[source].Acc[k];
if(SphP[place].MaxSignalVel < HydroDataPartialResult[source].MaxSignalVel)
SphP[place].MaxSignalVel = HydroDataPartialResult[source].MaxSignalVel;
}
}
}
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);
/* 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 update. A target
* particle is specified which may either be local, or reside in the
* communication buffer.
* Note that only the acceleration due to feedback and the change
* of signal velocity are updated, not the actual variation of Entropy rate.
*/
void get_sigvel_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, 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, vsig;
double h_j, dwk_j, r, r2, u, hfc_visc;
int phase=0;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
#ifdef ART_CONDUCTIVITY
printf("get_sigvel_evaluate is not implemented to run with the flag ART_CONDUCTIVITY");
endrun(674321);
#endif
#if defined(ART_VISCO_MM) || defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
printf("get_sigvel_evaluate is not implemented to run with the flags ART_VISCO_MM, ART_VISCO_RO, ART_VISCO_CD");
endrun(674322);
#endif
#ifdef DENSITY_INDEPENDENT_SPH
double egyrho, entvarpred;
#endif
if(mode == 0)
{
pos = P[target].Pos;
vel = SphP[target].VelPred;
h_i = SphP[target].Hsml;
mass = P[target].Mass;
dhsmlDensityFactor = SphP[target].DhsmlDensityFactor;
rho = SphP[target].Density;
#ifdef DENSITY_INDEPENDENT_SPH
pressure = updated_pressure(SphP[target].EntropyPred,SphP[target].EgyWtDensity,SphP[target].DeltaEgySpec);
#else
pressure = updated_pressure(SphP[target].EntropyPred,SphP[target].Density,SphP[target].DeltaEgySpec);
#endif
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);
#ifdef DENSITY_INDEPENDENT_SPH
egyrho = SphP[target].EgyWtDensity;
entvarpred = SphP[target].EntVarPred;
dhsmlDensityFactor = SphP[target].DhsmlEgyDensityFactor;
#endif
}
else
{
pos = HydroDataGet[target].Pos;
vel = HydroDataGet[target].Vel;
h_i = HydroDataGet[target].Hsml;
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;
#ifdef DENSITY_INDEPENDENT_SPH
egyrho = HydroDataGet[target].EgyRho;
entvarpred = HydroDataGet[target].EntVarPred;
#endif
}
/* initialize variables before SPH loop is started */
acc[0] = acc[1] = acc[2] = dtEntropy = 0;
maxSignalVel = 0;
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_i = pressure / (egyrho * egyrho);
#else
p_over_rho2_i = pressure / (rho * rho) * dhsmlDensityFactor;
#endif
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)
{
/*
here, we need to differenciate two cases, because DeltaEgySpec
is also used as a flag (see further) and may be equal to -1
*/
if (SphP[j].DeltaEgySpec>0) /* the particle is touch by feedback */
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_j = updated_pressure(SphP[j].EntropyPred,SphP[j].EgyWtDensity,SphP[j].DeltaEgySpec) / (SphP[j].EgyWtDensity * SphP[j].EgyWtDensity);
#else
p_over_rho2_j = updated_pressure(SphP[j].EntropyPred,SphP[j].Density,SphP[j].DeltaEgySpec) / (SphP[j].Density * SphP[j].Density);
#endif
else
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_j = SphP[j].Pressure / (SphP[j].EgyWtDensity * SphP[j].EgyWtDensity);
#else
p_over_rho2_j = SphP[j].Pressure / (SphP[j].Density * SphP[j].Density);
#endif
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_j = sqrt(GAMMA * SphP[j].Pressure / SphP[j].Density);
#else
soundspeed_j = sqrt(GAMMA * p_over_rho2_j * SphP[j].Density);
#endif
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;
if(P[j].Ti_endstep == All.Ti_Current)
if(vsig > SphP[j].MaxSignalVel)
SphP[j].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;
#ifndef DENSITY_INDEPENDENT_SPH
p_over_rho2_j *= SphP[j].DhsmlDensityFactor;
#endif
hfc_visc = 0.5 * P[j].Mass * visc * (dwk_i + dwk_j) / r;
#ifdef DENSITY_INDEPENDENT_SPH
hfc = hfc_visc;
/* leading-order term */
hfc += P[j].Mass *
(dwk_i*p_over_rho2_i*SphP[j].EntVarPred/entvarpred +
dwk_j*p_over_rho2_j*entvarpred/SphP[j].EntVarPred) / r;
/* grad-h corrections */
hfc += P[j].Mass *
(dwk_i*p_over_rho2_i*egyrho/rho*dhsmlDensityFactor +
dwk_j*p_over_rho2_j*SphP[j].EgyWtDensity/SphP[j].Density*SphP[j].DhsmlEgyDensityFactor) / r;
#else
hfc = hfc_visc + P[j].Mass * (p_over_rho2_i * dwk_i + p_over_rho2_j * dwk_j) / r;
#endif
acc[0] -= hfc * dx;
acc[1] -= hfc * dy;
acc[2] -= hfc * dz;
if(P[j].Ti_endstep == All.Ti_Current)
{
if(SphP[j].DeltaEgySpec == 0) /* the particle is active but not affected by feedback */
SphP[j].DeltaEgySpec = -1;
if(maxSignalVel > SphP[j].MaxSignalVel)
SphP[j].MaxSignalVel = maxSignalVel;
SphP[j].FeedbackUpdatedAccel[0] -= hfc * dx;
SphP[j].FeedbackUpdatedAccel[1] -= hfc * dy;
SphP[j].FeedbackUpdatedAccel[2] -= hfc * dz;
}
}
}
}
fflush(stdout);
}
while(startnode >= 0);
/* Now collect the result at the right place */
if(mode == 0)
{
SphP[target].MaxSignalVel = maxSignalVel;
for(k = 0; k < 3; k++)
SphP[target].FeedbackUpdatedAccel[k] = acc[k];
}
else
{
HydroDataResult[target].MaxSignalVel = maxSignalVel;
for(k = 0; k < 3; k++)
HydroDataResult[target].Acc[k] = acc[k];
}
}
/*! This function return the pressure as a function
* of entropy and density, but add the contribution
* of the feedback energy.
*/
FLOAT updated_pressure(FLOAT EntropyPred,FLOAT Density,FLOAT DeltaEgySpec)
{
FLOAT pressure;
FLOAT EgySpec,EgySpecUpdated;
/* energy from entropy */
EgySpec = EntropyPred / GAMMA_MINUS1 * pow(Density*a3inv, GAMMA_MINUS1);
EgySpecUpdated = EgySpec + DeltaEgySpec;
/* pressure */
pressure = GAMMA_MINUS1 * (Density*a3inv) * EgySpecUpdated;
return pressure;
}
/*! This function force the particle to be active
*/
void make_particle_active(int i)
{
int j;
int tstart, tend;
double dt_entr;
double dt_gravkick;
double dt_hydrokick;
//printf("(%d) make particle %d active.\n",ThisTask,i);
#ifdef PMGRID
double dt_gravkickA, dt_gravkickB;
#endif
tstart = (P[i].Ti_begstep + P[i].Ti_endstep) / 2; /* midpoint of old step */
tend = (P[i].Ti_begstep + All.Ti_Current) / 2; /* midpoint of shrinked step */
/* kick the particle back */
kickback(i,tstart,tend);
/* flag it as active */
P[i].Ti_endstep = All.Ti_Current;
NumForceUpdate++;
#ifdef PMGRID
printf("make_particle_active is not implemented to run with the flag PMGRID\n");
printf("here we need to add the last part of advance_and_find_timesteps()\n");
endrun(674323);
#endif
#ifdef MULTIPHASE
printf("make_particle_active is not implemented to run with the flag MULTIPHASE\n");
endrun(674324);
#endif
}
#endif
#endif
#ifdef SYNCHRONIZE_NGB_TIMESTEP
/*! In case we want to reduce the timestep of a particle,
* this function perform a negative kick on the particle.
* Positions, densities, predicted velocities and predicted
* entropies are not affected, as they are allready defined
* at the present step.
*/
void kickback(int i,int tstart,int tend)
{
int j;
double dt_entr;
double dt_gravkick;
double dt_hydrokick;
double DtEgySpec;
double a3inv;
#ifdef PMGRID
double dt_gravkickA, dt_gravkickB;
!!! this is probably bad, please check !!!
if(All.ComovingIntegrationOn)
dt_gravkickB = get_gravkick_factor(All.PM_Ti_begstep, All.Ti_Current) -
get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
else
dt_gravkickB = (All.Ti_Current - (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2) * All.Timebase_interval;
#endif
if(All.ComovingIntegrationOn)
{
dt_entr = (tend - tstart) * All.Timebase_interval;
dt_gravkick = get_gravkick_factor(tstart, tend);
dt_hydrokick = get_hydrokick_factor(tstart, tend);
a3inv = 1 / (All.Time * All.Time * All.Time);
}
else
{
dt_entr = dt_gravkick = dt_hydrokick = (tend - tstart) * All.Timebase_interval;
a3inv = 1;
}
for(j = 0; j < 3; j++)
{
P[i].Vel[j] += P[i].GravAccel[j] * dt_gravkick;
#ifdef PMGRID
P[i].Vel[j] += P[i].GravPM[j] * dt_gravkickB;
#endif
#ifdef AB_TURB
P[i].Vel[j] += SphP[i].TurbAccel[j] * dt_hydrokick;
#endif
#ifdef DISSIPATION_FORCES
P[i].Vel[j] += SphP[i].DissipationForcesAccel[j] * dt_hydrokick;
#endif
}
if(P[i].Type == 0) /* SPH stuff */
{
for(j = 0; j < 3; j++)
P[i].Vel[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
#ifdef COOLING
//DtEgySpec = - 1/GAMMA_MINUS1 * pow(SphP[i].Density * a3inv,GAMMA_MINUS1) * (SphP[i].DtEntropyRad);
DtEgySpec = SphP[i].DtEnergyRad; /* this is now independent of the density */
LocalSysState.RadiatedEnergy += DtEgySpec * dt_entr * P[i].Mass;
//printf("this is a check: dt_entr=%g RadiatedEnergy=%g\n",dt_entr,DtEgySpec * dt_entr * P[i].Mass);
#endif
SphP[i].Entropy += SphP[i].DtEntropy * dt_entr;
//#ifdef COOLING
// SphP[i].EntropyRad += -SphP[i].DtEntropyRad * dt_entr;
//#endif
}
#ifdef PMGRID
printf("kickback is not implemented to run with the flag PMGRID\n");
printf("here we need to add the last part of advance_and_find_timesteps()\n");
endrun(674324);
#endif
#ifdef MULTIPHASE
printf("kickback is not implemented to run with the flag MULTIPHASE\n");
endrun(674326);
#endif
#ifdef LIMIT_DVEL
printf("kickback is not implemented to run with the flag LIMIT_DVEL\n");
endrun(674327)
#endif
}
#endif /* SYNCHRONIZE_NGB_TIMESTEP */
c4science · Help