00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <string.h>
00004 #include <math.h>
00005 #include <mpi.h>
00006
00007 #include "allvars.h"
00008 #include "proto.h"
00009
00010
00022 void compute_global_quantities_of_system(void)
00023 {
00024 int i, j, n;
00025 struct state_of_system sys;
00026 double a1, a2, a3;
00027 double entr = 0, egyspec, vel[3];
00028 double dt_entr, dt_gravkick, dt_hydrokick;
00029
00030
00031
00032 if(All.ComovingIntegrationOn)
00033 {
00034 a1 = All.Time;
00035 a2 = All.Time * All.Time;
00036 a3 = All.Time * All.Time * All.Time;
00037 }
00038 else
00039 {
00040 a1 = a2 = a3 = 1;
00041 }
00042
00043
00044 for(n = 0; n < 6; n++)
00045 {
00046 sys.MassComp[n] = sys.EnergyKinComp[n] = sys.EnergyPotComp[n] = sys.EnergyIntComp[n] = 0;
00047
00048 for(j = 0; j < 4; j++)
00049 sys.CenterOfMassComp[n][j] = sys.MomentumComp[n][j] = sys.AngMomentumComp[n][j] = 0;
00050 }
00051
00052 for(i = 0; i < NumPart; i++)
00053 {
00054 sys.MassComp[P[i].Type] += P[i].Mass;
00055
00056 sys.EnergyPotComp[P[i].Type] += 0.5 * P[i].Mass * P[i].Potential / a1;
00057
00058 if(All.ComovingIntegrationOn)
00059 {
00060 dt_entr = (All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
00061 dt_gravkick = get_gravkick_factor(P[i].Ti_begstep, All.Ti_Current) -
00062 get_gravkick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
00063 dt_hydrokick = get_hydrokick_factor(P[i].Ti_begstep, All.Ti_Current) -
00064 get_hydrokick_factor(P[i].Ti_begstep, (P[i].Ti_begstep + P[i].Ti_endstep) / 2);
00065 }
00066 else
00067 dt_entr = dt_gravkick = dt_hydrokick =
00068 (All.Ti_Current - (P[i].Ti_begstep + P[i].Ti_endstep) / 2) * All.Timebase_interval;
00069
00070 for(j = 0; j < 3; j++)
00071 {
00072 vel[j] = P[i].Vel[j] + P[i].GravAccel[j] * dt_gravkick;
00073 if(P[i].Type == 0)
00074 vel[j] += SphP[i].HydroAccel[j] * dt_hydrokick;
00075 }
00076 if(P[i].Type == 0)
00077 entr = SphP[i].Entropy + SphP[i].DtEntropy * dt_entr;
00078
00079 #ifdef PMGRID
00080 if(All.ComovingIntegrationOn)
00081 dt_gravkick = get_gravkick_factor(All.PM_Ti_begstep, All.Ti_Current) -
00082 get_gravkick_factor(All.PM_Ti_begstep, (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2);
00083 else
00084 dt_gravkick = (All.Ti_Current - (All.PM_Ti_begstep + All.PM_Ti_endstep) / 2) * All.Timebase_interval;
00085
00086 for(j = 0; j < 3; j++)
00087 vel[j] += P[i].GravPM[j] * dt_gravkick;
00088 #endif
00089
00090 sys.EnergyKinComp[P[i].Type] +=
00091 0.5 * P[i].Mass * (vel[0] * vel[0] + vel[1] * vel[1] + vel[2] * vel[2]) / a2;
00092
00093 if(P[i].Type == 0)
00094 {
00095 #ifdef ISOTHERM_EQS
00096 egyspec = entr;
00097 #else
00098 egyspec = entr / (GAMMA_MINUS1) * pow(SphP[i].Density / a3, GAMMA_MINUS1);
00099 #endif
00100 sys.EnergyIntComp[0] += P[i].Mass * egyspec;
00101 }
00102
00103
00104
00105 for(j = 0; j < 3; j++)
00106 {
00107 sys.MomentumComp[P[i].Type][j] += P[i].Mass * vel[j];
00108 sys.CenterOfMassComp[P[i].Type][j] += P[i].Mass * P[i].Pos[j];
00109 }
00110
00111 sys.AngMomentumComp[P[i].Type][0] += P[i].Mass * (P[i].Pos[1] * vel[2] - P[i].Pos[2] * vel[1]);
00112 sys.AngMomentumComp[P[i].Type][1] += P[i].Mass * (P[i].Pos[2] * vel[0] - P[i].Pos[0] * vel[2]);
00113 sys.AngMomentumComp[P[i].Type][2] += P[i].Mass * (P[i].Pos[0] * vel[1] - P[i].Pos[1] * vel[0]);
00114 }
00115
00116
00117
00118 MPI_Reduce(&sys.MassComp[0], &SysState.MassComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00119 MPI_Reduce(&sys.EnergyPotComp[0], &SysState.EnergyPotComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00120 MPI_Reduce(&sys.EnergyIntComp[0], &SysState.EnergyIntComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00121 MPI_Reduce(&sys.EnergyKinComp[0], &SysState.EnergyKinComp[0], 6, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
00122 MPI_Reduce(&sys.MomentumComp[0][0], &SysState.MomentumComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
00123 MPI_COMM_WORLD);
00124 MPI_Reduce(&sys.AngMomentumComp[0][0], &SysState.AngMomentumComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
00125 MPI_COMM_WORLD);
00126 MPI_Reduce(&sys.CenterOfMassComp[0][0], &SysState.CenterOfMassComp[0][0], 6 * 4, MPI_DOUBLE, MPI_SUM, 0,
00127 MPI_COMM_WORLD);
00128
00129
00130 if(ThisTask == 0)
00131 {
00132 for(i = 0; i < 6; i++)
00133 SysState.EnergyTotComp[i] = SysState.EnergyKinComp[i] +
00134 SysState.EnergyPotComp[i] + SysState.EnergyIntComp[i];
00135
00136 SysState.Mass = SysState.EnergyKin = SysState.EnergyPot = SysState.EnergyInt = SysState.EnergyTot = 0;
00137
00138 for(j = 0; j < 3; j++)
00139 SysState.Momentum[j] = SysState.AngMomentum[j] = SysState.CenterOfMass[j] = 0;
00140
00141 for(i = 0; i < 6; i++)
00142 {
00143 SysState.Mass += SysState.MassComp[i];
00144 SysState.EnergyKin += SysState.EnergyKinComp[i];
00145 SysState.EnergyPot += SysState.EnergyPotComp[i];
00146 SysState.EnergyInt += SysState.EnergyIntComp[i];
00147 SysState.EnergyTot += SysState.EnergyTotComp[i];
00148
00149 for(j = 0; j < 3; j++)
00150 {
00151 SysState.Momentum[j] += SysState.MomentumComp[i][j];
00152 SysState.AngMomentum[j] += SysState.AngMomentumComp[i][j];
00153 SysState.CenterOfMass[j] += SysState.CenterOfMassComp[i][j];
00154 }
00155 }
00156
00157 for(i = 0; i < 6; i++)
00158 for(j = 0; j < 3; j++)
00159 if(SysState.MassComp[i] > 0)
00160 SysState.CenterOfMassComp[i][j] /= SysState.MassComp[i];
00161
00162 for(j = 0; j < 3; j++)
00163 if(SysState.Mass > 0)
00164 SysState.CenterOfMass[j] /= SysState.Mass;
00165
00166 for(i = 0; i < 6; i++)
00167 {
00168 SysState.CenterOfMassComp[i][3] = SysState.MomentumComp[i][3] = SysState.AngMomentumComp[i][3] = 0;
00169 for(j = 0; j < 3; j++)
00170 {
00171 SysState.CenterOfMassComp[i][3] +=
00172 SysState.CenterOfMassComp[i][j] * SysState.CenterOfMassComp[i][j];
00173 SysState.MomentumComp[i][3] += SysState.MomentumComp[i][j] * SysState.MomentumComp[i][j];
00174 SysState.AngMomentumComp[i][3] +=
00175 SysState.AngMomentumComp[i][j] * SysState.AngMomentumComp[i][j];
00176 }
00177 SysState.CenterOfMassComp[i][3] = sqrt(SysState.CenterOfMassComp[i][3]);
00178 SysState.MomentumComp[i][3] = sqrt(SysState.MomentumComp[i][3]);
00179 SysState.AngMomentumComp[i][3] = sqrt(SysState.AngMomentumComp[i][3]);
00180 }
00181
00182 SysState.CenterOfMass[3] = SysState.Momentum[3] = SysState.AngMomentum[3] = 0;
00183
00184 for(j = 0; j < 3; j++)
00185 {
00186 SysState.CenterOfMass[3] += SysState.CenterOfMass[j] * SysState.CenterOfMass[j];
00187 SysState.Momentum[3] += SysState.Momentum[j] * SysState.Momentum[j];
00188 SysState.AngMomentum[3] += SysState.AngMomentum[j] * SysState.AngMomentum[j];
00189 }
00190
00191 SysState.CenterOfMass[3] = sqrt(SysState.CenterOfMass[3]);
00192 SysState.Momentum[3] = sqrt(SysState.Momentum[3]);
00193 SysState.AngMomentum[3] = sqrt(SysState.AngMomentum[3]);
00194 }
00195
00196
00197 MPI_Bcast(&SysState, sizeof(struct state_of_system), MPI_BYTE, 0, MPI_COMM_WORLD);
00198 }
