Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F65311737
global.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
Sun, Jun 2, 20:30
Size
13 KB
Mime Type
text/x-c
Expires
Tue, Jun 4, 20:30 (2 d)
Engine
blob
Format
Raw Data
Handle
18002748
Attached To
rGEAR Gear
global.c
View Options
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include "allvars.h"
#include "proto.h"
/*! \file global.c
* \brief Computes global physical properties of the system
*/
/*! This routine computes various global properties of the particle
* distribution and stores the result in the struct `SysState'.
* Currently, not all the information that's computed here is actually
* used (e.g. momentum is not really used anywhere), just the energies are
* written to a log-file every once in a while.
*/
void compute_global_quantities_of_system(void)
{
int i, j, n;
struct state_of_system sys;
double a1, a2, a3;
double entr = 0, egyspec, vel[3];
#ifdef AGN_HEATING
double especagnheat = 0;
#endif
#ifdef SFR
/* see starformation.c */
#endif
#ifdef DISSIPATION_FORCES
double especdissipationforces = 0;
#endif
double dt_entr, dt_gravkick, dt_hydrokick;
if(All.ComovingIntegrationOn)
{
a1 = All.Time;
a2 = All.Time * All.Time;
a3 = All.Time * All.Time * All.Time;
}
else
{
a1 = a2 = a3 = 1;
}
for(n = 0; n < 6; n++)
{
sys.MassComp[n] = sys.EnergyKinComp[n] = sys.EnergyPotComp[n] = sys.EnergyIntComp[n] = 0;
#ifdef COOLING
sys.EnergyRadSphComp[n] = 0;
#endif
#ifdef AGN_HEATING
sys.EnergyAGNHeatComp[n] = 0;
#endif
#ifdef DISSIPATION_FORCES
sys.EnergyDissipationForcesComp[n] = 0;
#endif
for(j = 0; j < 4; j++)
sys.CenterOfMassComp[n][j] = sys.MomentumComp[n][j] = sys.AngMomentumComp[n][j] = 0;
}
#ifdef SFR
rearrange_particle_sequence();
#endif
for(i = 0; i < NumPart; i++)
{
sys.MassComp[P[i].Type] += P[i].Mass;
sys.EnergyPotComp[P[i].Type] += 0.5 * P[i].Mass * P[i].Potential / a1;
if(All.ComovingIntegrationOn)
{
dt_entr = (All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
dt_gravkick = get_gravkick_factor(P[i].Ti_begstep, All.Ti_Current) -
get_gravkick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
dt_hydrokick = get_hydrokick_factor(P[i].Ti_begstep, All.Ti_Current) -
get_hydrokick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
}
else
dt_entr = dt_gravkick = dt_hydrokick =
(All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
for(j = 0; j < 3; j++)
{
vel[j] = P[i].Vel[j] + P[i].GravAccel[j] * dt_gravkick;
if(P[i].Type == 0)
vel[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
#ifdef DISSIPATION_FORCES
if(P[i].Type == 0)
vel[j] += SphP[i].DissipationForcesAccel[j] * dt_hydrokick;
#endif
}
if(P[i].Type == 0)
{
entr = SphP[i].Entropy + SphP[i].DtEntropy * dt_entr;
#ifdef AGN_HEATING
especagnheat = SphP[i].EgySpecAGNHeat + SphP[i].DtEgySpecAGNHeat * dt_entr;
#endif
#ifdef DISSIPATION_FORCES
especdissipationforces = SphP[i].EnergyDissipationForces + SphP[i].DtEnergyDissipationForces * dt_entr;
#endif
#ifdef SFR
/* see starformation.c */
#endif
}
#ifdef PMGRID
if(All.ComovingIntegrationOn)
dt_gravkick = 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_gravkick = (All.Ti_Current - (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2) * All.Timebase_interval;
for(j = 0; j < 3; j++)
vel[j] += P[i].GravPM[j] * dt_gravkick;
#endif
sys.EnergyKinComp[P[i].Type] +=
0.5 * P[i].Mass * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]) / a2;
if(P[i].Type == 0)
{
if (entr>0)
#ifdef ISOTHERM_EQS
egyspec = entr;
#else
#ifdef MULTIPHASE
if (SphP[i].Phase== GAS_SPH)
#ifdef DENSITY_INDEPENDENT_SPH
egyspec = entr / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity / a3, GAMMA_MINUS1);
#else
egyspec = entr / (GAMMA_MINUS1) * pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
else
egyspec = entr;
#else
#ifdef DENSITY_INDEPENDENT_SPH
egyspec = entr / (GAMMA_MINUS1) * pow(SphP[i].EgyWtDensity / a3, GAMMA_MINUS1);
#else
egyspec = entr / (GAMMA_MINUS1) * pow(SphP[i].Density / a3, GAMMA_MINUS1);
#endif
#endif /*ISOTHERM_EQS*/
if(All.MinEgySpec)
egyspec = dmax(All.MinEgySpec,egyspec);
#endif
sys.EnergyIntComp[0] += P[i].Mass * egyspec;
#ifdef AGN_HEATING
sys.EnergyAGNHeatComp[0] += P[i].Mass * especagnheat;
#endif
#ifdef DISSIPATION_FORCES
sys.EnergyDissipationForcesComp[0] += P[i].Mass * especdissipationforces;
#endif
#ifdef SFR
/* see starformation.c */
#endif
}
for(j = 0; j < 3; j++)
{
sys.MomentumComp[P[i].Type][j] += P[i].Mass * vel[j];
sys.CenterOfMassComp[P[i].Type][j] += P[i].Mass * P[i].Pos[j];
}
sys.AngMomentumComp[P[i].Type][0] += P[i].Mass * (P[i].Pos[1] * vel[2] - P[i].Pos[2] * vel[1]);
sys.AngMomentumComp[P[i].Type][1] += P[i].Mass * (P[i].Pos[2] * vel[0] - P[i].Pos[0] * vel[2]);
sys.AngMomentumComp[P[i].Type][2] += P[i].Mass * (P[i].Pos[0] * vel[1] - P[i].Pos[1] * vel[0]);
}
/* count energy lost by different processes */
#ifdef COOLING
sys.EnergyRadSphComp[0]=LocalSysState.RadiatedEnergy;
#endif
#ifdef SFR
sys.EnergyIntComp[ST]=LocalSysState.StarEnergyInt;
#endif
/* count thermal feedback */
#ifdef CHIMIE_THERMAL_FEEDBACK
sys.EnergyThermalFeedbackComp[0]=LocalSysState.EnergyThermalFeedback;
for(i = 1; i < 6; i++)
sys.EnergyThermalFeedbackComp[i]=0;
#endif
/* count kinetic feedback */
#ifdef CHIMIE_KINETIC_FEEDBACK
sys.EnergyKineticFeedbackComp[0]=LocalSysState.EnergyKineticFeedback;
for(i = 1; i < 6; i++)
sys.EnergyKineticFeedbackComp[i]=0;
#endif
/* count sticky */
#ifdef MULTIPHASE
sys.EnergyRadStickyComp[0]=LocalSysState.EnergyRadSticky;
for(i = 1; i < 6; i++)
sys.EnergyRadStickyComp[i]=0;
#endif
/* count feedback wind */
#ifdef FEEDBACK_WIND
sys.EnergyFeedbackWindComp[0]=LocalSysState.EnergyFeedbackWind;
for(i = 1; i < 6; i++)
sys.EnergyFeedbackWindComp[i]=0;
#endif
/* count dissipated energy */
#ifdef INTEGRAL_CONSERVING_DISSIPATION
sys.EnergyICDissipationComp[0]=LocalSysState.EnergyICDissipation;
for(i = 1; i < 6; i++)
sys.EnergyICDissipationComp[i]=0;
#endif
/* some the stuff over all processors */
MPI_Reduce(&sys.MassComp[0], &SysState.MassComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&sys.EnergyPotComp[0], &SysState.EnergyPotComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&sys.EnergyIntComp[0], &SysState.EnergyIntComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#ifdef COOLING
MPI_Reduce(&sys.EnergyRadSphComp[0], &SysState.EnergyRadSphComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef AGN_HEATING
MPI_Reduce(&sys.EnergyAGNHeatComp[0], &SysState.EnergyAGNHeatComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef DISSIPATION_FORCES
MPI_Reduce(&sys.EnergyDissipationForcesComp[0], &SysState.EnergyDissipationForcesComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
MPI_Reduce(&sys.EnergyThermalFeedbackComp[0], &SysState.EnergyThermalFeedbackComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
MPI_Reduce(&sys.EnergyKineticFeedbackComp[0], &SysState.EnergyKineticFeedbackComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef MULTIPHASE
MPI_Reduce(&sys.EnergyRadStickyComp[0], &SysState.EnergyRadStickyComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef FEEDBACK_WIND
MPI_Reduce(&sys.EnergyFeedbackWindComp[0], &SysState.EnergyFeedbackWindComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
MPI_Reduce(&sys.EnergyICDissipationComp[0], &SysState.EnergyICDissipationComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
#endif
MPI_Reduce(&sys.EnergyKinComp[0], &SysState.EnergyKinComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&sys.MomentumComp[0][0], &SysState.MomentumComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&sys.AngMomentumComp[0][0], &SysState.AngMomentumComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD);
MPI_Reduce(&sys.CenterOfMassComp[0][0], &SysState.CenterOfMassComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
MPI_COMM_WORLD);
#ifdef BUBBLES
SysState.EnergyBubblesComp[0]=All.EnergyBubbles;
for(i = 1; i < 6; i++)
SysState.EnergyBubblesComp[i]=0;
#endif
if(ThisTask == 0)
{
for(i = 0; i < 6; i++)
{
SysState.EnergyTotComp[i] = SysState.EnergyKinComp[i]
+ SysState.EnergyPotComp[i]
+ SysState.EnergyIntComp[i];
#ifdef COOLING
SysState.EnergyTotComp[i] += SysState.EnergyRadSphComp[i];
#endif
#ifdef AGN_HEATING
SysState.EnergyTotComp[i] += SysState.EnergyAGNHeatComp[i];
#endif
#ifdef DISSIPATION_FORCES
SysState.EnergyTotComp[i] += SysState.EnergyDissipationForcesComp[i];
#endif
#ifdef MULTIPHASE
SysState.EnergyTotComp[i] += SysState.EnergyRadStickyComp[i];
#endif
#ifdef FEEDBACK_WIND
SysState.EnergyTotComp[i] += SysState.EnergyFeedbackWindComp[i];
#endif
#ifdef BUBBLES
SysState.EnergyTotComp[i] += SysState.EnergyBubblesComp[i];
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
SysState.EnergyTotComp[i] += SysState.EnergyThermalFeedbackComp[i];
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
SysState.EnergyTotComp[i] += SysState.EnergyKineticFeedbackComp[i];
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
SysState.EnergyTotComp[i] += SysState.EnergyICDissipationComp[i];
#endif
}
SysState.Mass = SysState.EnergyKin = SysState.EnergyPot = SysState.EnergyInt = SysState.EnergyTot = 0;
#ifdef COOLING
SysState.EnergyRadSph = 0;
#endif
#ifdef AGN_HEATING
SysState.EnergyAGNHeat = 0;
#endif
#ifdef DISSIPATION_FORCES
SysState.EnergyDissipationForces = 0;
#endif
#ifdef MULTIPHASE
SysState.EnergyRadSticky = 0;
#endif
#ifdef FEEDBACK_WIND
SysState.EnergyFeedbackWind = 0;
#endif
#ifdef BUBBLES
SysState.EnergyBubbles = 0;
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
SysState.EnergyThermalFeedback = 0;
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
SysState.EnergyKineticFeedback = 0;
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
SysState.EnergyICDissipation = 0;
#endif
for(j = 0; j < 3; j++)
SysState.Momentum[j] = SysState.AngMomentum[j] = SysState.CenterOfMass[j] = 0;
for(i = 0; i < 6; i++)
{
SysState.Mass += SysState.MassComp[i];
SysState.EnergyKin += SysState.EnergyKinComp[i];
SysState.EnergyPot += SysState.EnergyPotComp[i];
SysState.EnergyInt += SysState.EnergyIntComp[i];
SysState.EnergyTot += SysState.EnergyTotComp[i];
#ifdef COOLING
SysState.EnergyRadSph += SysState.EnergyRadSphComp[i];
#endif
#ifdef AGN_HEATING
SysState.EnergyAGNHeat += SysState.EnergyAGNHeatComp[i];
#endif
#ifdef DISSIPATION_FORCES
SysState.EnergyDissipationForces += SysState.EnergyDissipationForcesComp[i];
#endif
#ifdef MULTIPHASE
SysState.EnergyRadSticky += SysState.EnergyRadStickyComp[i];
#endif
#ifdef FEEDBACK_WIND
SysState.EnergyFeedbackWind += SysState.EnergyFeedbackWindComp[i];
#endif
#ifdef BUBBLES
SysState.EnergyBubbles += SysState.EnergyBubblesComp[i];
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
SysState.EnergyThermalFeedback += SysState.EnergyThermalFeedbackComp[i];
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
SysState.EnergyKineticFeedback += SysState.EnergyKineticFeedbackComp[i];
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
SysState.EnergyICDissipation += SysState.EnergyICDissipationComp[i];
#endif
for(j = 0; j < 3; j++)
{
SysState.Momentum[j] += SysState.MomentumComp[i][j];
SysState.AngMomentum[j] += SysState.AngMomentumComp[i][j];
SysState.CenterOfMass[j] += SysState.CenterOfMassComp[i][j];
}
}
for(i = 0; i < 6; i++)
for(j = 0; j < 3; j++)
if(SysState.MassComp[i] > 0)
SysState.CenterOfMassComp[i][j] /= SysState.MassComp[i];
for(j = 0; j < 3; j++)
if(SysState.Mass > 0)
SysState.CenterOfMass[j] /= SysState.Mass;
for(i = 0; i < 6; i++)
{
SysState.CenterOfMassComp[i][3] = SysState.MomentumComp[i][3] = SysState.AngMomentumComp[i][3] = 0;
for(j = 0; j < 3; j++)
{
SysState.CenterOfMassComp[i][3] +=
SysState.CenterOfMassComp[i][j] * SysState.CenterOfMassComp[i][j];
SysState.MomentumComp[i][3] += SysState.MomentumComp[i][j] * SysState.MomentumComp[i][j];
SysState.AngMomentumComp[i][3] +=
SysState.AngMomentumComp[i][j] * SysState.AngMomentumComp[i][j];
}
SysState.CenterOfMassComp[i][3] = sqrt(SysState.CenterOfMassComp[i][3]);
SysState.MomentumComp[i][3] = sqrt(SysState.MomentumComp[i][3]);
SysState.AngMomentumComp[i][3] = sqrt(SysState.AngMomentumComp[i][3]);
}
SysState.CenterOfMass[3] = SysState.Momentum[3] = SysState.AngMomentum[3] = 0;
for(j = 0; j < 3; j++)
{
SysState.CenterOfMass[3] += SysState.CenterOfMass[j] * SysState.CenterOfMass[j];
SysState.Momentum[3] += SysState.Momentum[j] * SysState.Momentum[j];
SysState.AngMomentum[3] += SysState.AngMomentum[j] * SysState.AngMomentum[j];
}
SysState.CenterOfMass[3] = sqrt(SysState.CenterOfMass[3]);
SysState.Momentum[3] = sqrt(SysState.Momentum[3]);
SysState.AngMomentum[3] = sqrt(SysState.AngMomentum[3]);
}
/* give everyone the result, maybe the want to do something with it */
MPI_Bcast(&SysState, sizeof(struct state_of_system), MPI_BYTE, 0, MPI_COMM_WORLD);
}
Event Timeline
Log In to Comment