Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65173878
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
Sat, Jun 1, 12:27
Size
56 KB
Mime Type
text/x-c
Expires
Mon, Jun 3, 12:27 (2 d)
Engine
blob
Format
Raw Data
Handle
18018644
Attached To
rGEAR Gear
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 ART_VISCO_CD
int ii,jj;
#endif
#ifdef DETAILED_CPU_OUTPUT_IN_HYDRA
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
#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
}
#ifdef OUTPUT_CONDUCTIVITY
for (i=0;i<N_gas;i++)
SphP[i].OptVar1 = sqrt( pow(SphP[i].GradEnergyInt[0],2)+pow(SphP[i].GradEnergyInt[1],2)+pow(SphP[i].GradEnergyInt[2],2))*SphP[i].Hsml/(SphP[i].Pressure/(SphP[i].Density*GAMMA_MINUS1)) ;
#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++;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK)
for(j = 0; j < 3; j++)
SphP[n].FeedbackUpdatedAccel[j] = 0;
#endif
#ifdef METAL_DIFFUSION
SphP[n].MetalDiffusionA = 0;
for(j = 0; j < NELEMENTS; j++)
SphP[n].MetalDiffusionB[j] = 0;
#endif
}
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;
#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;
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
#ifdef ART_CONDUCTIVITY
HydroDataIn[nexport].NormGradEnergyInt = sqrt( pow(SphP[i].GradEnergyInt[0],2)+pow(SphP[i].GradEnergyInt[1],2)+pow(SphP[i].GradEnergyInt[2],2));
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
HydroDataIn[nexport].ArtBulkViscConst = SphP[i].ArtBulkViscConst;
#endif
/* calculation of F1 */
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_i = sqrt(GAMMA * SphP[i].Pressure / SphP[i].EgyWtDensity);
#else
soundspeed_i = sqrt(GAMMA * SphP[i].Pressure / SphP[i].Density);
#endif
HydroDataIn[nexport].F1 = fabs(SphP[i].DivVel) /
(fabs(SphP[i].DivVel) + SphP[i].CurlVel +
0.0001 * soundspeed_i / SphP[i].Hsml / fac_mu);
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[i].DeltaEgySpec>0) /* the particle is touch by feedback */
{
#ifdef DENSITY_INDEPENDENT_SPH
HydroDataIn[nexport].PressureFeedbackUpdated = updated_pressure_hydra(SphP[i].EntropyPred,SphP[i].EgyWtDensity,SphP[i].DeltaEgySpec);
#else
HydroDataIn[nexport].PressureFeedbackUpdated = updated_pressure_hydra(SphP[i].EntropyPred,SphP[i].Density,SphP[i].DeltaEgySpec);
#endif
}
else
{
HydroDataIn[nexport].PressureFeedbackUpdated = SphP[i].Pressure;
}
/* calculation of F1 */
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_i = sqrt(GAMMA * HydroDataIn[nexport].PressureFeedbackUpdated / SphP[i].EgyWtDensity);
#else
soundspeed_i = sqrt(GAMMA * HydroDataIn[nexport].PressureFeedbackUpdated / SphP[i].Density);
#endif
HydroDataIn[nexport].F1FeedbackUpdated = fabs(SphP[i].DivVel) /
(fabs(SphP[i].DivVel) + SphP[i].CurlVel +
0.0001 * soundspeed_i / SphP[i].Hsml / fac_mu);
#endif
#ifdef METAL_DIFFUSION
HydroDataIn[nexport].MetalDiffusionCoeff = SphP[i].MetalDiffusionCoeff;
#endif
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
#ifdef OUTPUT_CONDUCTIVITY
SphP[place].OptVar2 += HydroDataPartialResult[source].OptVar2;
#endif
#ifdef ART_VISCO_CD
/* reduce DmatCD and TmatCD*/
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
SphP[place].DmatCD[ii][jj] += HydroDataPartialResult[source].DmatCD[ii][jj];
SphP[place].TmatCD[ii][jj] += HydroDataPartialResult[source].TmatCD[ii][jj];
}
SphP[place].R_CD += HydroDataPartialResult[source].R_CD;
if(SphP[place].MaxSignalVelCD < HydroDataPartialResult[source].MaxSignalVelCD)
SphP[place].MaxSignalVelCD = HydroDataPartialResult[source].MaxSignalVelCD;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
for(k = 0; k < 3; k++)
SphP[place].FeedbackUpdatedAccel[k] += HydroDataPartialResult[source].AccFeedbackUpdated[k];
if(SphP[place].MaxSignalVelFeedbackUpdated < HydroDataPartialResult[source].maxSignalVelFeedbackUpdated)
SphP[place].MaxSignalVelFeedbackUpdated = HydroDataPartialResult[source].maxSignalVelFeedbackUpdated;
#endif
#ifdef METAL_DIFFUSION
for(k = 0; k < NELEMENTS; k++)
SphP[place].MetalDiffusionB[k] += HydroDataPartialResult[source].MetalDiffusionB[k];
SphP[place].MetalDiffusionA += HydroDataPartialResult[source].MetalDiffusionA;
#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();
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
for(i = 0; i < N_gas; i++)
SphP[i].MaxSignalVel=dmax(SphP[i].MaxSignalVel,SphP[i].MaxSignalVelFeedbackUpdated);
#endif
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
{
#ifdef DENSITY_INDEPENDENT_SPH
SphP[i].DtEntropy *= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].EgyWtDensity, GAMMA_MINUS1));
#else
SphP[i].DtEntropy *= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].Density, GAMMA_MINUS1));
#endif
#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].OptVar1 = sqrt(v2m - v1m*v1m);
else
SphP[i].OptVar1 = 0.0;
}
else
SphP[i].OptVar1 = 0.0;
#endif
#ifdef OUTPUT_CONDUCTIVITY
SphP[i].OptVar2*= GAMMA_MINUS1 / (hubble_a2 * pow(SphP[i].Density, GAMMA_MINUS1)); /* to dA/dt */
if (SphP[i].OptVar2!=0)
SphP[i].OptVar2=SphP[i].Entropy/fabs(SphP[i].OptVar2); /* to time*/
else
SphP[i].OptVar2=0;
#endif
#ifdef ART_VISCO_CD
compute_art_visc(i);
#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;
}
#ifdef DETAILED_CPU_OUTPUT_IN_HYDRA
numlist = malloc(sizeof(int) * NTask);
timecomplist = malloc(sizeof(double) * NTask);
timecommsummlist = malloc(sizeof(double) * NTask);
timeimbalancelist = malloc(sizeof(double) * NTask);
MPI_Gather(&NumSphUpdate, 1, MPI_INT, numlist, 1, MPI_INT, 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 hydra\n\n");
fprintf(FdTimings, "Nupdate ");
for (i=0;i<NTask;i++)
fprintf(FdTimings, "%12d ",numlist[i]); /* nombre de part par proc */
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(numlist);
#endif
}
/*! 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 WITH_ID_IN_HYDRA
int id;
#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
#ifdef ART_CONDUCTIVITY
double hfc_cond,vsig_u,vsig_u_max;
double Arho_i,Arho_j;
double u_i,u_j;
double normGradEnergyInt_i, normGradEnergyInt_j;
double dtEntropy_artcond;
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)|| defined(ART_VISCO_CD)
double alpha_i, alpha_j;
double alpha_ij;
double beta_ij;
double soundspeed_ij;
#ifdef ART_VISCO_CD
double wk,wk_i,wk_j,hinv3;
int ii, jj;
double DmatCD[3][3];
double TmatCD[3][3];
double dv_dot_dx;
double R_CD;
double sign_DiVel;
double maxSignalVelCD;
#endif
#endif
#ifdef DENSITY_INDEPENDENT_SPH
double egyrho, entvarpred;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
double p_over_rho2_iFeedbackUpdated,soundspeed_iFeedbackUpdated;
double p_over_rho2_jFeedbackUpdated,soundspeed_jFeedbackUpdated;
double maxSignalVelFeedbackUpdated,hfcFeedbackUpdated,hfc_viscFeedbackUpdated,viscFeedbackUpdated;
double accFeedbackUpdated[3];
FLOAT f1FeedbackUpdated,f2FeedbackUpdated;
FLOAT pressureFeedbackUpdated;
#ifdef DENSITY_INDEPENDENT_SPH
double pup;
#endif
#endif
#ifdef METAL_DIFFUSION
FLOAT MetalDiffusionCoeff;
double Kij=0;
double MetalDiffusionA=0;
double MetalDiffusionB[NELEMENTS];
for (k = 0;k < NELEMENTS; k++)
MetalDiffusionB[k]=0;
#endif
if(mode == 0)
{
pos = P[target].Pos;
vel = SphP[target].VelPred;
h_i = SphP[target].Hsml;
#ifdef FEEDBACK
EnergySN = SphP[target].EnergySN;
#endif
#ifdef WITH_ID_IN_HYDRA
id = P[target].ID;
#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;
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_i = sqrt(GAMMA * pressure / SphP[target].EgyWtDensity);
#else
soundspeed_i = sqrt(GAMMA * pressure / rho);
#endif
f1 = fabs(SphP[target].DivVel) /
(fabs(SphP[target].DivVel) + SphP[target].CurlVel +
0.0001 * soundspeed_i / SphP[target].Hsml / fac_mu);
#ifdef ART_CONDUCTIVITY
normGradEnergyInt_i = sqrt( pow(SphP[target].GradEnergyInt[0],2)+pow(SphP[target].GradEnergyInt[1],2)+pow(SphP[target].GradEnergyInt[2],2));
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO) || defined(ART_VISCO_CD)
alpha_i = SphP[target].ArtBulkViscConst;
#endif
#ifdef DENSITY_INDEPENDENT_SPH
egyrho = SphP[target].EgyWtDensity;
entvarpred = SphP[target].EntVarPred;
dhsmlDensityFactor = SphP[target].DhsmlEgyDensityFactor;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[target].DeltaEgySpec>0)
{
#ifdef DENSITY_INDEPENDENT_SPH
pressureFeedbackUpdated = updated_pressure_hydra(SphP[target].EntropyPred,SphP[target].EgyWtDensity,SphP[target].DeltaEgySpec);
#else
pressureFeedbackUpdated = updated_pressure_hydra(SphP[target].EntropyPred,SphP[target].Density,SphP[target].DeltaEgySpec);
#endif
}
else
{
pressureFeedbackUpdated = SphP[target].Pressure;
}
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_iFeedbackUpdated = sqrt(GAMMA * pressureFeedbackUpdated / SphP[target].EgyWtDensity);
#else
soundspeed_iFeedbackUpdated = sqrt(GAMMA * pressureFeedbackUpdated / rho);
#endif
f1FeedbackUpdated = fabs(SphP[target].DivVel) /
(fabs(SphP[target].DivVel) + SphP[target].CurlVel +
0.0001 * soundspeed_iFeedbackUpdated / SphP[target].Hsml / fac_mu);
#endif
#ifdef METAL_DIFFUSION
MetalDiffusionCoeff = SphP[target].MetalDiffusionCoeff;
#endif
}
else
{
pos = HydroDataGet[target].Pos;
vel = HydroDataGet[target].Vel;
h_i = HydroDataGet[target].Hsml;
#ifdef FEEDBACK
EnergySN = HydroDataGet[target].EnergySN;
#endif
#ifdef WITH_ID_IN_HYDRA
id = HydroDataGet[target].ID;
#endif
mass = HydroDataGet[target].Mass;
dhsmlDensityFactor = HydroDataGet[target].DhsmlDensityFactor;
rho = HydroDataGet[target].Density;
pressure = HydroDataGet[target].Pressure;
timestep = HydroDataGet[target].Timestep;
#ifdef DENSITY_INDEPENDENT_SPH
soundspeed_i = sqrt(GAMMA * pressure / HydroDataGet[target].EgyRho);
#else
soundspeed_i = sqrt(GAMMA * pressure / rho);
#endif
f1 = HydroDataGet[target].F1;
#ifdef ART_CONDUCTIVITY
normGradEnergyInt_i = HydroDataGet[target].NormGradEnergyInt;
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)|| defined(ART_VISCO_CD)
alpha_i = HydroDataGet[target].ArtBulkViscConst;
#endif
#ifdef DENSITY_INDEPENDENT_SPH
egyrho = HydroDataGet[target].EgyRho;
entvarpred = HydroDataGet[target].EntVarPred;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
pressureFeedbackUpdated = HydroDataGet[target].PressureFeedbackUpdated;
#endif
#ifdef METAL_DIFFUSION
MetalDiffusionCoeff = HydroDataGet[target].MetalDiffusionCoeff;
#endif
}
/* initialize variables before SPH loop is started */
acc[0] = acc[1] = acc[2] = dtEntropy = 0;
maxSignalVel = 0;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
accFeedbackUpdated[0] = accFeedbackUpdated[1] = accFeedbackUpdated[2] = 0;
maxSignalVelFeedbackUpdated = 0;
#endif
#ifdef FEEDBACK
dtEgySpecFeedback=0;
#endif
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_i = pressure / (egyrho * egyrho);
#else
p_over_rho2_i = pressure / (rho * rho) * dhsmlDensityFactor;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_iFeedbackUpdated = pressureFeedbackUpdated / (egyrho * egyrho);
#else
p_over_rho2_iFeedbackUpdated = pressureFeedbackUpdated / (rho * rho) * dhsmlDensityFactor;
#endif
#endif
h_i2 = h_i * h_i;
#ifdef ART_CONDUCTIVITY
Arho_i = pressure/rho;
u_i = pressure/(rho*GAMMA_MINUS1);
dtEntropy_artcond=0;
#endif
#ifdef ART_VISCO_CD
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
DmatCD[ii][jj] = 0.;
TmatCD[ii][jj] = 0.;
}
R_CD = 0.;
maxSignalVelCD=0;
#endif
/* 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)
{
#ifdef DENSITY_INDEPENDENT_SPH
p_over_rho2_j = SphP[j].Pressure / (SphP[j].EgyWtDensity * SphP[j].EgyWtDensity);
soundspeed_j = sqrt(GAMMA * SphP[j].Pressure / SphP[j].EgyWtDensity);
#else
p_over_rho2_j = SphP[j].Pressure / (SphP[j].Density * SphP[j].Density);
soundspeed_j = sqrt(GAMMA * p_over_rho2_j * SphP[j].Density);
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[j].DeltaEgySpec>0) /* the particle is touch by feedback */
{
#ifdef DENSITY_INDEPENDENT_SPH
pup = updated_pressure_hydra(SphP[j].EntropyPred,SphP[j].EgyWtDensity,SphP[j].DeltaEgySpec);
p_over_rho2_jFeedbackUpdated = pup / (SphP[j].EgyWtDensity * SphP[j].EgyWtDensity);
soundspeed_jFeedbackUpdated = sqrt(GAMMA * pup / SphP[j].EgyWtDensity);
#else
p_over_rho2_jFeedbackUpdated = updated_pressure_hydra(SphP[j].EntropyPred,SphP[j].Density,SphP[j].DeltaEgySpec) / (SphP[j].Density * SphP[j].Density);
soundspeed_jFeedbackUpdated = sqrt(GAMMA * p_over_rho2_jFeedbackUpdated * SphP[j].Density);
#endif
}
else
{
p_over_rho2_jFeedbackUpdated = p_over_rho2_j;
soundspeed_jFeedbackUpdated = soundspeed_j;
}
#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 defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[j].DeltaEgySpec>0) /* the particle is touch by feedback */
{
if(soundspeed_iFeedbackUpdated + soundspeed_jFeedbackUpdated > maxSignalVelFeedbackUpdated)
maxSignalVelFeedbackUpdated = soundspeed_iFeedbackUpdated + soundspeed_jFeedbackUpdated;
}
#endif
/*********************************/
/* standard form of viscosity */
/*********************************/
#if !defined(ART_VISCO_MM) && !defined(ART_VISCO_RO) && !defined(ART_VISCO_CD)
visc = artificial_viscosity(r,vdotr2,soundspeed_i,soundspeed_j,dwk_i,dwk_j,timestep,j,mass,rho,f1,&maxSignalVel);
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[j].DeltaEgySpec>0)
viscFeedbackUpdated = artificial_viscosity(r,vdotr2,soundspeed_iFeedbackUpdated,soundspeed_jFeedbackUpdated,dwk_i,dwk_j,timestep,j,mass,rho,f1FeedbackUpdated,&maxSignalVelFeedbackUpdated);
#endif
#endif
/**************************************/
/* alternative form of viscosity RO MM*/
/**************************************/
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)
visc = artificial_viscosity_improved(r,vdotr2,soundspeed_i,soundspeed_j,dwk_i,dwk_j,h_i,alpha_i,timestep,j,mass,rho,f1,&maxSignalVel);
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[j].DeltaEgySpec>0)
viscFeedbackUpdated = artificial_viscosity_improved(r,vdotr2,soundspeed_iFeedbackUpdated,soundspeed_jFeedbackUpdated,dwk_i,dwk_j,h_i,alpha_i,timestep,j,mass,rho,f1FeedbackUpdated,&maxSignalVelFeedbackUpdated);
#endif
#endif
/**************************************/
/* alternative form of viscosity CD */
/**************************************/
#if defined(ART_VISCO_CD)
visc = artificial_viscosity_CD(r,vdotr2,soundspeed_i,soundspeed_j,dwk_i,dwk_j,timestep,j,mass,rho,f1,&maxSignalVel);
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (SphP[j].DeltaEgySpec>0)
viscFeedbackUpdated = artificial_viscosity_CD(r,vdotr2,soundspeed_iFeedbackUpdated,soundspeed_jFeedbackUpdated,dwk_i,dwk_j,timestep,j,mass,rho,f1FeedbackUpdated,&maxSignalVelFeedbackUpdated);
#endif
#endif
/*********************************/
/* now finish sph */
/*********************************/
#ifndef DENSITY_INDEPENDENT_SPH
p_over_rho2_j *= SphP[j].DhsmlDensityFactor;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
p_over_rho2_jFeedbackUpdated *= SphP[j].DhsmlDensityFactor;
#endif
#endif
hfc_visc = 0.5 * P[j].Mass * visc * (dwk_i + dwk_j) / r;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
hfc_viscFeedbackUpdated = 0.5 * P[j].Mass * viscFeedbackUpdated * (dwk_i + dwk_j) / r;
#endif
#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;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
hfcFeedbackUpdated = hfc_viscFeedbackUpdated;
/* leading-order term */
hfcFeedbackUpdated += P[j].Mass *
(dwk_i*p_over_rho2_iFeedbackUpdated*SphP[j].EntVarPred/entvarpred +
dwk_j*p_over_rho2_jFeedbackUpdated*entvarpred/SphP[j].EntVarPred) / r;
/* grad-h corrections */
hfcFeedbackUpdated += P[j].Mass *
(dwk_i*p_over_rho2_iFeedbackUpdated*egyrho/rho*dhsmlDensityFactor +
dwk_j*p_over_rho2_jFeedbackUpdated*SphP[j].EgyWtDensity/SphP[j].Density*SphP[j].DhsmlEgyDensityFactor) / r;
#endif
#else /* if not DENSITY_INDEPENDENT_SPH */
hfc = hfc_visc + P[j].Mass * (p_over_rho2_i * dwk_i + p_over_rho2_j * dwk_j) / r;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
hfcFeedbackUpdated = hfc_viscFeedbackUpdated + P[j].Mass * (p_over_rho2_iFeedbackUpdated * dwk_i + p_over_rho2_jFeedbackUpdated * dwk_j) / r;
#endif
#endif /* DENSITY_INDEPENDENT_SPH */
acc[0] -= hfc * dx;
acc[1] -= hfc * dy;
acc[2] -= hfc * dz;
dtEntropy += 0.5 * hfc_visc * vdotr2;
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
accFeedbackUpdated[0] -= hfcFeedbackUpdated * dx;
accFeedbackUpdated[1] -= hfcFeedbackUpdated * dy;
accFeedbackUpdated[2] -= hfcFeedbackUpdated * 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(maxSignalVelFeedbackUpdated > SphP[j].MaxSignalVelFeedbackUpdated)
{
SphP[j].MaxSignalVelFeedbackUpdated = maxSignalVelFeedbackUpdated;
}
SphP[j].FeedbackUpdatedAccel[0] -= hfcFeedbackUpdated * dx;
SphP[j].FeedbackUpdatedAccel[1] -= hfcFeedbackUpdated * dy;
SphP[j].FeedbackUpdatedAccel[2] -= hfcFeedbackUpdated * dz;
}
#endif
/*********************************/
/* prediction for the next step */
/*********************************/
#ifdef ART_VISCO_CD
artificial_viscosity_CD_prediction(r,vdotr2,soundspeed_i,soundspeed_j,dwk_i,dwk_j,timestep,j,mass,rho,f1,&maxSignalVel);
#endif
/*****************************************/
/* ART CONDUCTIVITY */
/*****************************************/
#ifdef ART_CONDUCTIVITY
mu_ij = fac_mu * vdotr2 / r;
vsig_u = soundspeed_i + soundspeed_j - 3 * mu_ij;
Arho_j = SphP[j].Pressure / SphP[j].Density;
rho_ij = 0.5 * (rho + SphP[j].Density);
/* switch */
normGradEnergyInt_j = sqrt( pow(SphP[j].GradEnergyInt[0],2)+pow(SphP[j].GradEnergyInt[1],2)+pow(SphP[j].GradEnergyInt[2],2));
u_j = SphP[j].Pressure/(SphP[j].Density*GAMMA_MINUS1);
hfc_cond = P[j].Mass * All.ArtCondConst * vsig_u * (Arho_i-Arho_j)/ rho_ij * 0.5*(dwk_i + dwk_j);
dtEntropy_artcond = hfc_cond / GAMMA_MINUS1;
dtEntropy += dtEntropy_artcond ;
#endif
/*****************************************/
/* METAL DIFFUSION */
/*****************************************/
#ifdef METAL_DIFFUSION
if(MetalDiffusionCoeff!=0 && SphP[j].MetalDiffusionCoeff!=0)
{
Kij = P[j].Mass / (SphP[j].Density * rho ) * 4.0 * MetalDiffusionCoeff * SphP[j].MetalDiffusionCoeff / ( MetalDiffusionCoeff + SphP[j].MetalDiffusionCoeff ) * 0.5 * (dwk_i + dwk_j) / r ;
MetalDiffusionA += Kij;
for(k=0;k<NELEMENTS;k++)
MetalDiffusionB[k] += Kij * SphP[j].Metal[k] ;
}
#endif
/*****************************************/
/* 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
#ifdef OUTPUT_CONDUCTIVITY
SphP[target].OptVar2 = dtEntropy_artcond;
#endif
#ifdef ART_VISCO_CD
/*collect DmatCD, TmatCD*/
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
SphP[target].DmatCD[ii][jj] = DmatCD[ii][jj];
SphP[target].TmatCD[ii][jj] = TmatCD[ii][jj];
}
SphP[target].R_CD = R_CD;
SphP[target].MaxSignalVelCD = maxSignalVelCD;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
if (maxSignalVelFeedbackUpdated>SphP[target].MaxSignalVelFeedbackUpdated)
SphP[target].MaxSignalVelFeedbackUpdated = maxSignalVelFeedbackUpdated;
for(k = 0; k < 3; k++)
SphP[target].FeedbackUpdatedAccel[k] = accFeedbackUpdated[k];
#endif
#ifdef METAL_DIFFUSION
SphP[target].MetalDiffusionA = MetalDiffusionA;
for(k = 0; k < NELEMENTS; k++)
SphP[target].MetalDiffusionB[k] = MetalDiffusionB[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
#ifdef OUTPUT_CONDUCTIVITY
HydroDataResult[target].OptVar2 = dtEntropy_artcond;
#endif
#ifdef ART_VISCO_CD
/*collect DmatCD, TmatCD*/
for (ii = 0; ii < 3; ii++)
for (jj = 0; jj < 3; jj++)
{
HydroDataResult[target].DmatCD[ii][jj] = DmatCD[ii][jj];
HydroDataResult[target].TmatCD[ii][jj] = TmatCD[ii][jj];
}
HydroDataResult[target].R_CD = R_CD;
HydroDataResult[target].MaxSignalVelCD = maxSignalVelCD;
#endif
#if defined(TIMESTEP_UPDATE_FOR_FEEDBACK) && defined(CHIMIE_THERMAL_FEEDBACK)
HydroDataResult[target].maxSignalVelFeedbackUpdated = maxSignalVelFeedbackUpdated;
for(k = 0; k < 3; k++)
HydroDataResult[target].AccFeedbackUpdated[k] = accFeedbackUpdated[k];
#endif
#ifdef METAL_DIFFUSION
HydroDataResult[target].MetalDiffusionA = MetalDiffusionA;
for(k = 0; k < NELEMENTS; k++)
HydroDataResult[target].MetalDiffusionB[k] = MetalDiffusionB[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;
}
/*! default from Gadget-2
*/
double artificial_viscosity(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
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;
return visc;
}
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)
/*! Tiret version
*/
double artificial_viscosity_improved_old_tiret(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
double alpha_i,alpha_j,alpha_ij,beta_ij,soundspeed_ij;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
if(vdotr2 < 0) /* ... artificial viscosity */
{
alpha_j = SphP[j].ArtBulkViscConst;
alpha_ij = 0.5*(alpha_i + alpha_j);
beta_ij = 3/2. * alpha_ij; /* 3/2 is compatible with Springel 05 */
mu_ij = fac_mu * vdotr2 / r; /* note: this is negative! */
vsig = soundspeed_i + soundspeed_j - 2*beta_ij/alpha_ij * mu_ij;
if(vsig > *maxSignalVel)
*maxSignalVel = vsig;
soundspeed_ij = 0.5 * (soundspeed_i + soundspeed_j);
f2 =
fabs(SphP[j].DivVel) / (fabs(SphP[j].DivVel) + SphP[j].CurlVel +
0.0001 * soundspeed_j / fac_mu / SphP[j].Hsml);
rho_ij = 0.5 * (rho + SphP[j].Density);
visc = (- alpha_ij * soundspeed_ij * mu_ij + beta_ij * mu_ij * mu_ij) / rho_ij * 0.5*(f1 + f2) ;
#ifndef NOVISCOSITYLIMITER
if(vdotr2 < 0)
{
/* 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;
return visc;
}
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)
/*! Monaghan + variation of coefficient (!!! this may diverge for small rij)
*/
double artificial_viscosity_improved_old_Arnaudon(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,double h_i,double alpha_i,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
double alpha_j,alpha_ij,beta_ij,soundspeed_ij;
double h_ij,r2;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
if(vdotr2 < 0) /* ... artificial viscosity */
{
r2 = r*r;
h_ij= 0.5 * (h_i + SphP[j].Hsml);
mu_ij = h_ij* fac_mu * vdotr2 / r2;
alpha_j = SphP[j].ArtBulkViscConst;
alpha_ij = 0.5*(alpha_i + alpha_j);
beta_ij = 2 * alpha_ij; /* 3/2 is compatible with Springel 05 */
soundspeed_ij = 0.5 * (soundspeed_i + soundspeed_j);
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 = (-alpha_ij *soundspeed_ij + beta_ij*mu_ij )*mu_ij / rho_ij *0.5 *(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;
return visc;
}
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)
/*! Rosswog (with divergence corrected)
*/
double artificial_viscosity_improved(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,double h_i,double alpha_i,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
double alpha_j,alpha_ij,beta_ij,soundspeed_ij;
double h_ij,r2;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
if(vdotr2 < 0) /* ... artificial viscosity */
{
r2 = r*r;
h_ij= 0.5 * (h_i + SphP[j].Hsml);
mu_ij = h_ij* fac_mu * vdotr2 / (r2+0.01*(h_ij*h_ij));
vsig = soundspeed_i + soundspeed_j - 3 * mu_ij;
if(vsig > *maxSignalVel)
*maxSignalVel = vsig;
alpha_j = SphP[j].ArtBulkViscConst;
alpha_ij = 0.5*(alpha_i + alpha_j);
beta_ij = 2 * alpha_ij; /* 3/2 is compatible with Springel 05 */
soundspeed_ij = 0.5 * (soundspeed_i + soundspeed_j);
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 = (-alpha_ij *soundspeed_ij + beta_ij*mu_ij )*mu_ij / rho_ij *0.5 *(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;
return visc;
}
#endif
#if defined(ART_VISCO_MM)|| defined(ART_VISCO_RO)
/*! Springel + Rosswog
*/
double artificial_viscosity_improved_springel(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,double h_i,double alpha_i,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
double alpha_j,alpha_ij,beta_ij,soundspeed_ij;
double h_ij,r2;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
if(vdotr2 < 0) /* ... artificial viscosity */
{
alpha_j = SphP[j].ArtBulkViscConst;
alpha_ij = 0.5*(alpha_i + alpha_j);
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 * alpha_ij * 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;
return visc;
}
#endif
#if defined(ART_VISCO_CD)
/*! Cullen-Dehnen artificial viscosity
*/
double artificial_viscosity_CD(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
double visc,vsig,rho_ij,mu_ij;
FLOAT f2;
#ifndef NOVISCOSITYLIMITER
double dt;
#endif
alpha_j = SphP[j].ArtBulkViscConst;
alpha_ij = 0.5*(alpha_i + alpha_j);
beta_ij = 3/2. * alpha_ij;
mu_ij = fac_mu * vdotr2 / r; /* note: this is negative! */
vsig = soundspeed_i + soundspeed_j - 2*beta_ij/alpha_ij * mu_ij;
if(vsig > maxSignalVel)
maxSignalVel = vsig;
soundspeed_ij = 0.5 * (soundspeed_i + soundspeed_j);
rho_ij = 0.5 * (rho + SphP[j].Density);
visc = (- alpha_ij * soundspeed_ij * mu_ij + beta_ij * mu_ij * mu_ij) / rho_ij;
#ifndef NOVISCOSITYLIMITER
if(vdotr2 < 0)
{
/* 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
return visc;
}
#endif
#if defined(ART_VISCO_CD)
/*! Cullen-Dehnen artificial viscosity
*/
double artificial_viscosity_CD_prediction(double r,double vdotr2,double soundspeed_i,double soundspeed_j,double dwk_i,double dwk_j,double mu_ij,int timestep,int j,FLOAT mass,FLOAT rho,FLOAT f1,double *maxSignalVel)
{
/* WARNING WARNING WARNING WARNING WARNING WARNING WARNING
This part is not finished and should not compile.
WARNING WARNING WARNING WARNING WARNING WARNING WARNING * /
/* COMPUTE wk_i, wk_j, wk */
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); */
wk = 0.5*(dwk_i + dwk_j)/r;
/* CHOICE OF THE WEIGHT */
wk = P[j].Mass * wk / SphP[j].Density;
/* COMPUTE the matrix Di, Ti */
DmatCD[0][0] += dvx * dx * wk;
DmatCD[1][0] += dvy * dx * wk;
DmatCD[2][0] += dvz * dx * wk;
DmatCD[0][1] += dvx * dy * wk;
DmatCD[1][1] += dvy * dy * wk;
DmatCD[2][1] += dvz * dy * wk;
DmatCD[0][2] += dvx * dz * wk;
DmatCD[1][2] += dvy * dz * wk;
DmatCD[2][2] += dvz * dz * wk;
TmatCD[0][0] += dx * dx * wk;
TmatCD[1][0] += dy * dx * wk;
TmatCD[2][0] += dz * dx * wk;
TmatCD[0][1] += dx * dy * wk;
TmatCD[1][1] += dy * dy * wk;
TmatCD[2][1] += dz * dy * wk;
TmatCD[0][2] += dx * dz * wk;
TmatCD[1][2] += dy * dz * wk;
TmatCD[2][2] += dz * dz * wk;
/* COMPUTE maxSignalVel */
dv_dot_dx = dvx * dx + dvy * dy + dvz * dz;
vsig = soundspeed_ij - dmin(0., dv_dot_dx);
if(vsig > maxSignalVelCD)
maxSignalVelCD = vsig;
/* compute chock indicator */
if (SphP[j].DiVelAccurate>0.)
sign_DiVel = 1.;
else
sign_DiVel = -1.;
R_CD += 0.5 * sign_DiVel * P[j].Mass * (wk_i + wk_j);
}
#endif
#ifdef TIMESTEP_UPDATE_FOR_FEEDBACK
/*! This function return the pressure as a function
* of entropy and density, but add the contribution
* of the feedback energy.
*/
FLOAT updated_pressure_hydra(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;
}
#endif
#ifdef JEANS_PRESSURE_FLOOR
double JeansPressureFloor(int i)
{
double Pjeans;
double density;
#ifdef DENSITY_INDEPENDENT_SPH
density = SphP[i].EgyWtDensity;
#else
density = SphP[i].Density;
#endif
Pjeans = 4/PI * All.JeansMassFactor * (SphP[i].Hsml*SphP[i].Hsml) * (density*density) * All.G/GAMMA;
return Pjeans;
}
#endif
Event Timeline
Log In to Comment