Homec4science
Diffusion Gear (master) Browse

Gear/src/run.c
b87f34a0ce7emaster

Last Change
Raw File

run.c
View Options

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <mpi.h>
#include <unistd.h>
#include "allvars.h"
#include "proto.h"
/*! \file run.c
* \brief iterates over timesteps, main loop
*/
/*! This routine contains the main simulation loop that iterates over single
* timesteps. The loop terminates when the cpu-time limit is reached, when a
* `stop' file is found in the output directory, or when the simulation ends
* because we arrived at TimeMax.
*/
void run(void)
{
FILE *fd;
int stopflag = 0;
char stopfname[200], contfname[200];
double t0, t1;
#ifdef DETAILED_CPU
double tstart,tend;
#endif
sprintf(stopfname, "%sstop", All.OutputDir);
sprintf(contfname, "%scont", All.OutputDir);
unlink(contfname);
do /* main loop */
{
t0 = second();
find_next_sync_point_and_drift(); /* find next synchronization point and drift particles to this time.
* If needed, this function will also write an output file
* at the desired time.
*/
every_timestep_stuff(); /* write some info to log-files */
#ifdef PNBODY
compute_pnbody();
#endif
#ifdef OUTPUT_EVERY_TIMESTEP
savepositions(All.SnapshotFileCount++); /* write snapshot file */
#endif
#ifdef DETAILED_CPU
tstart = second();
#endif
#ifdef AGN_ACCRETION
compute_agn_accretion(); /* compute accretion */
#endif
#ifdef BONDI_ACCRETION
compute_bondi_accretion(); /* compute bondi accretion */
#endif
#ifdef BUBBLES
make_bubble(); /* create a bubble */
#endif
#ifdef MULTIPHASE
update_phase(); /* allow particles to change their phase */
#endif
#ifdef CORIOLIS
set_outer_potential_coriolis(); /* coriolis */
#endif
#ifdef CHIMIE
chimie();
#endif
#ifdef COOLING_FCT_FROM_HDF5
checkRedshiftForUpdate();
#endif
#ifdef COOLING_WIERSMA
float a = get_a_from_CosmicTime(All.Time);
float Redshift = get_Redshift_from_a(a);
setTablesFromRedshift(Redshift);
#endif
#ifdef FOF
fof();
#else // FOF
#ifdef SFR
#ifdef COMPUTE_SFR_ENERGY
density(1);
force_update_hmax();
sfr_compute_energy_int(1);
#endif
star_formation(); /* starformation */
#endif // SFR
#endif // FOF
#ifdef GAS_ACCRETION
gas_accretion();
#endif
#ifdef HOT_HALO
halo_creation();
#endif
//#ifdef VANISHING_PARTICLES
// vanishing_particles();
//#endif
#ifdef MULTIPHASE
sticky();
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
integral_conserving_dissipation();
#endif
#ifdef DETAILED_CPU
tend = second();
All.CPU_Physics += timediff(tstart, tend);
#endif
domain_Decomposition(); /* do domain decomposition if needed */
compute_accelerations(0); /* compute accelerations for
* the particles that are to be advanced
*/
#if defined(SFR) && defined(COMPUTE_SFR_ENERGY)
#ifdef DETAILED_CPU
tstart = second();
#endif
//sfr_compute_energy_int(2);
#ifdef DETAILED_CPU
tend = second();
All.CPU_Physics += timediff(tstart, tend);
#endif
#endif
/* check whether we want a full energy statistics */
if((All.Time - All.TimeLastStatistics) >= All.TimeBetStatistics)
{
#ifdef COMPUTE_POTENTIAL_ENERGY
compute_potential();
#endif
#ifndef ADVANCEDSTATISTICS
energy_statistics(); /* compute and output energy statistics */
#else
advanced_energy_statistics(); /* compute and output energy statistics */
#endif
All.TimeLastStatistics += All.TimeBetStatistics;
}
advance_and_find_timesteps(); /* 'kick' active particles in
* momentum space and compute new
* timesteps for them
*/
#ifdef SIDM_ON
scatter(); /* Find neighbours, calculate scatter rates
* and kick those particles that scatter */
#endif
#ifdef PERIODICOUTER
move_tracer(All.TimeStep);
#ifdef FLEXSTEPS
kick_tracer(All.PresentMinStep*All.Timebase_interval/2. + All.TimeStep/2.);
update_tracequat(All.PresentMinStep*All.Timebase_interval/2. + All.TimeStep/2.);
#else
if(All.Time - All.TimeStep == All.TimeBegin)
{
kick_tracer(All.TimeStep/2.);
update_tracequat(All.TimeStep/2.);
}
else
{
kick_tracer(All.TimeStep);
update_tracequat(All.TimeStep);
}
#endif
#endif
All.NumCurrentTiStep++;
/* Check whether we need to interrupt the run */
if(ThisTask == 0)
{
/* Is the stop-file present? If yes, interrupt the run. */
if((fd = fopen(stopfname, "r")))
{
fclose(fd);
stopflag = 1;
unlink(stopfname);
}
/* are we running out of CPU-time ? If yes, interrupt run. */
//if(CPUThisRun > 0.85 * All.TimeLimitCPU)
if(CPUThisRun > 1.0 * All.TimeLimitCPU)
{
printf("reaching time-limit. stopping.\n");
stopflag = 2;
}
}
MPI_Bcast(&stopflag, 1, MPI_INT, 0, MPI_COMM_WORLD);
if(stopflag)
{
restart(0); /* write restart file */
MPI_Barrier(MPI_COMM_WORLD);
if(stopflag == 2 && ThisTask == 0)
{
if((fd = fopen(contfname, "w")))
fclose(fd);
}
if(stopflag == 2 && All.ResubmitOn && ThisTask == 0)
{
close_outputfiles();
system(All.ResubmitCommand);
}
return;
}
/* is it time to write a regular restart-file? (for security) */
if(ThisTask == 0)
{
if((CPUThisRun - All.TimeLastRestartFile) >= All.CpuTimeBetRestartFile)
{
All.TimeLastRestartFile = CPUThisRun;
stopflag = 3;
}
else
stopflag = 0;
}
MPI_Bcast(&stopflag, 1, MPI_INT, 0, MPI_COMM_WORLD);
if(stopflag == 3)
{
restart(0); /* write an occasional restart file */
stopflag = 0;
}
t1 = second();
All.CPU_Total += timediff(t0, t1);
CPUThisRun += timediff(t0, t1);
}
while(All.Ti_Current < TIMEBASE && All.Time <= All.TimeMax);
restart(0);
savepositions(All.SnapshotFileCount++); /* write a last snapshot
* file at final time (will
* be overwritten if
* All.TimeMax is increased
* and the run is continued)
*/
}
/*! This function finds the next synchronization point of the system (i.e. the
* earliest point of time any of the particles needs a force computation),
* and drifts the system to this point of time. If the system drifts over
* the desired time of a snapshot file, the function will drift to this
* moment, generate an output, and then resume the drift.
*/
void find_next_sync_point_and_drift(void)
{
int n, min, min_glob, flag, *temp;
double timeold;
double t0, t1;
#ifdef DETAILED_CPU
double tstart,tend;
#endif
t0 = second();
#ifdef DETAILED_CPU
tstart = t0;
#endif
timeold = All.Time;
for(n = 1, min = P[0].Ti_endstep; n < NumPart; n++)
if(min > P[n].Ti_endstep)
min = P[n].Ti_endstep;
MPI_Allreduce(&min, &min_glob, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
/* We check whether this is a full step where all particles are synchronized */
flag = 1;
for(n = 0; n < NumPart; n++)
if(P[n].Ti_endstep > min_glob)
flag = 0;
MPI_Allreduce(&flag, &Flag_FullStep, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
#ifdef PMGRID
if(min_glob >= All.PM_Ti_endstep)
{
min_glob = All.PM_Ti_endstep;
Flag_FullStep = 1;
}
#endif
/* Determine 'NumForceUpdate', i.e. the number of particles on this processor that are going to be active */
for(n = 0, NumForceUpdate = 0; n < NumPart; n++)
{
if(P[n].Ti_endstep == min_glob)
#ifdef SELECTIVE_NO_GRAVITY
if(!((1 << P[n].Type) & (SELECTIVE_NO_GRAVITY)))
#endif
NumForceUpdate++;
}
/* note: NumForcesSinceLastDomainDecomp has type "long long" */
temp = malloc(NTask * sizeof(int));
MPI_Allgather(&NumForceUpdate, 1, MPI_INT, temp, 1, MPI_INT, MPI_COMM_WORLD);
for(n = 0; n < NTask; n++)
All.NumForcesSinceLastDomainDecomp += temp[n];
#ifdef COUNT_ACTIVE_PARTICLES
long long NumActivePatricles;
NumActivePatricles=0;
for(n = 0; n < NTask; n++)
NumActivePatricles+=temp[n];
#endif
free(temp);
t1 = second();
All.CPU_Predict += timediff(t0, t1);
while(min_glob >= All.Ti_nextoutput && All.Ti_nextoutput >= 0)
{
move_particles(All.Ti_Current, All.Ti_nextoutput);
All.Ti_Current = All.Ti_nextoutput;
if(All.ComovingIntegrationOn)
All.Time = All.TimeBegin * exp(All.Ti_Current * All.Timebase_interval);
else
All.Time = All.TimeBegin + All.Ti_Current * All.Timebase_interval;
#ifdef OUTPUTPOTENTIAL
All.NumForcesSinceLastDomainDecomp = 1 + All.TotNumPart * All.TreeDomainUpdateFrequency;
domain_Decomposition();
compute_potential();
#endif
#ifndef OUTPUT_EVERY_TIMESTEP
savepositions(All.SnapshotFileCount++); /* write snapshot file */
#endif
All.Ti_nextoutput = find_next_outputtime(All.Ti_nextoutput + 1);
}
move_particles(All.Ti_Current, min_glob);
All.Ti_Current = min_glob;
if(All.ComovingIntegrationOn)
All.Time = All.TimeBegin * exp(All.Ti_Current * All.Timebase_interval);
else
All.Time = All.TimeBegin + All.Ti_Current * All.Timebase_interval;
All.TimeStep = All.Time - timeold;
#ifdef COUNT_ACTIVE_PARTICLES
if (ThisTask==0)
{
fprintf(FdTimings,"===========================================================================================\n");
fprintf(FdTimings,"Step = %06d Time = %g \n\n",All.NumCurrentTiStep,All.Time);
fprintf(FdTimings,"%g %g : Total number of active particles : %d%09d\n\n",All.Time,All.TimeStep,(int) (NumActivePatricles / 1000000000), (int) (NumActivePatricles % 1000000000));
}
int *numpartlist;
int i;
int tot;
numpartlist = malloc(sizeof(int) * NTask*6);
MPI_Gather(&N_gas, 1, MPI_INT, &numpartlist[NTask*0], 1, MPI_INT, 0, MPI_COMM_WORLD);
#ifdef STELLAR_PROP
MPI_Gather(&N_stars, 1, MPI_INT, &numpartlist[NTask*1], 1, MPI_INT, 0, MPI_COMM_WORLD);
#endif
MPI_Gather(&NumPart, 1, MPI_INT, &numpartlist[NTask*2], 1, MPI_INT, 0, MPI_COMM_WORLD);
if (ThisTask==0)
{
tot = 0;
fprintf(FdTimings,"gas ");
for (i=0;i<NTask;i++)
{
fprintf(FdTimings, "%12d ",numpartlist[NTask*0+i]); /* nombre de part par proc */
tot += numpartlist[NTask*0+i];
}
fprintf(FdTimings, " : %12d ",tot);
fprintf(FdTimings,"\n");
tot = 0;
fprintf(FdTimings,"stars ");
for (i=0;i<NTask;i++)
{
fprintf(FdTimings, "%12d ",numpartlist[NTask*1+i]); /* nombre de part par proc */
tot += numpartlist[NTask*1+i];
}
fprintf(FdTimings, " : %12d ",tot);
fprintf(FdTimings,"\n");
tot = 0;
fprintf(FdTimings,"remaining ");
for (i=0;i<NTask;i++)
{
fprintf(FdTimings, "%12d ",numpartlist[NTask*2+i]-numpartlist[NTask*1+i]-numpartlist[NTask*0+i]); /* nombre de part par proc */
tot += numpartlist[NTask*2+i]-numpartlist[NTask*1+i]-numpartlist[NTask*0+i];
}
fprintf(FdTimings, " : %12d ",tot);
fprintf(FdTimings,"\n\n");
tot = 0;
fprintf(FdTimings,"total ");
for (i=0;i<NTask;i++)
{
fprintf(FdTimings, "%12d ",numpartlist[NTask*2+i]); /* nombre de part par proc */
tot += numpartlist[NTask*2+i];
}
fprintf(FdTimings, " : %12d ",tot);
fprintf(FdTimings,"\n\n");
fflush(FdTimings);
}
free(numpartlist);
#endif
#ifdef DETAILED_CPU
tend = second();
All.CPU_Leapfrog += timediff(tstart, tend);
#endif
}
/*! this function returns the next output time that is equal or larger to
* ti_curr
*/
int find_next_outputtime(int ti_curr)
{
int i, ti, ti_next, iter = 0;
double next, time;
ti_next = -1;
if(All.OutputListOn)
{
for(i = 0; i < All.OutputListLength; i++)
{
time = All.OutputListTimes[i];
if(time >= All.TimeBegin && time <= All.TimeMax)
{
if(All.ComovingIntegrationOn)
ti = log(time / All.TimeBegin) / All.Timebase_interval;
else
ti = (time - All.TimeBegin) / All.Timebase_interval;
if(ti >= ti_curr)
{
if(ti_next == -1)
ti_next = ti;
if(ti_next > ti)
ti_next = ti;
}
}
}
}
else
{
if(All.ComovingIntegrationOn)
{
if(All.TimeBetSnapshot <= 1.0)
{
printf("TimeBetSnapshot > 1.0 required for your simulation.\n");
endrun(13123);
}
}
else
{
if(All.TimeBetSnapshot <= 0.0)
{
printf("TimeBetSnapshot > 0.0 required for your simulation.\n");
endrun(13123);
}
}
time = All.TimeOfFirstSnapshot;
iter = 0;
while(time < All.TimeBegin)
{
if(All.ComovingIntegrationOn)
time *= All.TimeBetSnapshot;
else
time += All.TimeBetSnapshot;
iter++;
if(iter > 1000000)
{
printf("Can't determine next output time.\n");
endrun(110);
}
}
while(time <= All.TimeMax)
{
if(All.ComovingIntegrationOn)
ti = log(time / All.TimeBegin) / All.Timebase_interval;
else
ti = (time - All.TimeBegin) / All.Timebase_interval;
if(ti >= ti_curr)
{
ti_next = ti;
break;
}
if(All.ComovingIntegrationOn)
time *= All.TimeBetSnapshot;
else
time += All.TimeBetSnapshot;
iter++;
if(iter > 1000000)
{
printf("Can't determine next output time.\n");
endrun(111);
}
}
}
if(ti_next == -1)
{
ti_next = 2 * TIMEBASE; /* this will prevent any further output */
if(ThisTask == 0)
printf("\nThere is no valid time for a further snapshot file.\n");
}
else
{
if(All.ComovingIntegrationOn)
next = All.TimeBegin * exp(ti_next * All.Timebase_interval);
else
next = All.TimeBegin + ti_next * All.Timebase_interval;
if(ThisTask == 0)
printf("\nSetting next time for snapshot file to Time_next= %g\n\n", next);
}
return ti_next;
}
/*! This routine writes one line for every timestep to two log-files. In
* FdInfo, we just list the timesteps that have been done, while in FdCPU the
* cumulative cpu-time consumption in various parts of the code is stored.
*/
void every_timestep_stuff(void)
{
double z;
#ifdef DETAILED_CPU
double tstart,tend;
tstart = second();
#endif
if(ThisTask == 0)
{
if(All.ComovingIntegrationOn)
{
z = 1.0 / (All.Time) - 1;
fprintf(FdInfo, "\nBegin Step %d, Time: %g, Redshift: %g, Systemstep: %g, Dloga: %g\n",
All.NumCurrentTiStep, All.Time, z, All.TimeStep,
log(All.Time) - log(All.Time - All.TimeStep));
printf("\nBegin Step %d, Time: %g, Redshift: %g, Systemstep: %g, Dloga: %g\n", All.NumCurrentTiStep,
All.Time, z, All.TimeStep, log(All.Time) - log(All.Time - All.TimeStep));
fflush(FdInfo);
}
else
{
fprintf(FdInfo, "\nBegin Step %d, Time: %g, Systemstep: %g\n", All.NumCurrentTiStep, All.Time,
All.TimeStep);
printf("\nBegin Step %d, Time: %g, Systemstep: %g\n", All.NumCurrentTiStep, All.Time, All.TimeStep);
fflush(FdInfo);
}
printf("-------------------------------------------------------------\n");
fflush(stdout);
#ifdef ADVANCEDCPUSTATISTICS
fprintf(FdCPU, "%d ", All.NumCurrentTiStep);
fprintf(FdCPU, "%g ", All.Time);
fprintf(FdCPU, "%d ", NTask);
fprintf(FdCPU,"%10.2f ",All.CPU_Total);
#ifdef DETAILED_CPU
fprintf(FdCPU,"%10.2f ",All.CPU_Leapfrog);
fprintf(FdCPU,"%10.2f ",All.CPU_Physics);
fprintf(FdCPU,"%10.2f ",All.CPU_Residual);
fprintf(FdCPU,"%10.2f ",All.CPU_Accel);
fprintf(FdCPU,"%10.2f ",All.CPU_Begrun);
#endif
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity);
fprintf(FdCPU,"%10.2f ",All.CPU_Hydro);
#ifdef COOLING
fprintf(FdCPU,"%10.2f ",All.CPU_Cooling);
#endif
#ifdef SFR
fprintf(FdCPU,"%10.2f ",All.CPU_StarFormation);
#endif
#ifdef CHIMIE
fprintf(FdCPU,"%10.2f ",All.CPU_Chimie);
#endif
#ifdef MULTIPHASE
fprintf(FdCPU,"%10.2f ",All.CPU_Sticky);
#endif
fprintf(FdCPU,"%10.2f ",All.CPU_Domain);
fprintf(FdCPU,"%10.2f ",All.CPU_Potential);
fprintf(FdCPU,"%10.2f ",All.CPU_Predict);
fprintf(FdCPU,"%10.2f ",All.CPU_TimeLine);
fprintf(FdCPU,"%10.2f ",All.CPU_Snapshot);
fprintf(FdCPU,"%10.2f ",All.CPU_TreeWalk);
fprintf(FdCPU,"%10.2f ",All.CPU_TreeConstruction);
fprintf(FdCPU,"%10.2f ",All.CPU_CommSum);
fprintf(FdCPU,"%10.2f ",All.CPU_Imbalance);
fprintf(FdCPU,"%10.2f ",All.CPU_HydCompWalk);
fprintf(FdCPU,"%10.2f ",All.CPU_HydCommSumm);
fprintf(FdCPU,"%10.2f ",All.CPU_HydImbalance);
fprintf(FdCPU,"%10.2f ",All.CPU_EnsureNgb);
fprintf(FdCPU,"%10.2f ",All.CPU_PM);
fprintf(FdCPU,"%10.2f ",All.CPU_Peano);
#ifdef DETAILED_CPU_DOMAIN
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_findExtend);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_determineTopTree);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_sumCost);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_findSplit);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_shiftSplit);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_countToGo);
fprintf(FdCPU,"%10.2f ",All.CPU_Domain_exchange);
#endif
#ifdef DETAILED_CPU_GRAVITY
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_TreeWalk1);
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_TreeWalk2);
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_CommSum1);
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_CommSum2);
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_Imbalance1);
fprintf(FdCPU,"%10.2f ",All.CPU_Gravity_Imbalance2);
#endif
fprintf(FdCPU,"\n");
fflush(FdCPU);
#else
fprintf(FdCPU, "Step %d, Time: %g, CPUs: %d\n", All.NumCurrentTiStep, All.Time, NTask);
fprintf(FdCPU,
"%10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f\n",
All.CPU_Total, All.CPU_Gravity, All.CPU_Hydro, All.CPU_Domain, All.CPU_Potential,
All.CPU_Predict, All.CPU_TimeLine, All.CPU_Snapshot, All.CPU_TreeWalk, All.CPU_TreeConstruction,
All.CPU_CommSum, All.CPU_Imbalance, All.CPU_HydCompWalk, All.CPU_HydCommSumm,
All.CPU_HydImbalance, All.CPU_EnsureNgb, All.CPU_PM, All.CPU_Peano);
fflush(FdCPU);
#endif
}
set_random_numbers();
#ifdef DETAILED_CPU
tend = second();
All.CPU_Residual += timediff(tstart, tend);
#endif
}
/*! This routine first calls a computation of various global quantities of the
* particle distribution, and then writes some statistics about the energies
* in the various particle components to the file FdEnergy.
*/
void energy_statistics(void)
{
compute_global_quantities_of_system();
if(ThisTask == 0)
{
fprintf(FdEnergy,
"%g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g %g\n",
All.Time, SysState.EnergyInt, SysState.EnergyPot, SysState.EnergyKin, SysState.EnergyIntComp[0],
SysState.EnergyPotComp[0], SysState.EnergyKinComp[0], SysState.EnergyIntComp[1],
SysState.EnergyPotComp[1], SysState.EnergyKinComp[1], SysState.EnergyIntComp[2],
SysState.EnergyPotComp[2], SysState.EnergyKinComp[2], SysState.EnergyIntComp[3],
SysState.EnergyPotComp[3], SysState.EnergyKinComp[3], SysState.EnergyIntComp[4],
SysState.EnergyPotComp[4], SysState.EnergyKinComp[4], SysState.EnergyIntComp[5],
SysState.EnergyPotComp[5], SysState.EnergyKinComp[5], SysState.MassComp[0],
SysState.MassComp[1], SysState.MassComp[2], SysState.MassComp[3], SysState.MassComp[4],
SysState.MassComp[5]);
fflush(FdEnergy);
}
}
/*! This routine first calls a computation of various global quantities of the
* particle distribution, and then writes some statistics about the energies
* in the various particle components to the file FdEnergy.
*/
#ifdef ADVANCEDSTATISTICS
void advanced_energy_statistics(void)
{
int i;
#ifdef DETAILED_CPU
double tstart,tend;
tstart = second();
#endif
compute_global_quantities_of_system();
if(ThisTask == 0)
{
/**************/
/* energy */
/**************/
/* time */
fprintf(FdEnergy,"%g ",All.Time);
/* total */
fprintf(FdEnergy,"%g %g %g ",SysState.EnergyInt, SysState.EnergyPot, SysState.EnergyKin);
#ifdef COOLING
fprintf(FdEnergy,"%g ",SysState.EnergyRadSph);
#endif
#ifdef AGN_HEATING
fprintf(FdEnergy,"%g ",SysState.EnergyAGNHeat);
#endif
#ifdef DISSIPATION_FORCES
fprintf(FdEnergy,"%g ",SysState.EnergyDissipationForces);
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
fprintf(FdEnergy,"%g ",SysState.EnergyICDissipation);
#endif
#ifdef MULTIPHASE
fprintf(FdEnergy,"%g ",SysState.EnergyRadSticky);
#endif
#ifdef FEEDBACK_WIND
fprintf(FdEnergy,"%g ",SysState.EnergyFeedbackWind);
#endif
#ifdef BUBBLES
fprintf(FdEnergy,"%g ",SysState.EnergyBubbles);
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
fprintf(FdEnergy,"%g ",SysState.EnergyThermalFeedback);
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
fprintf(FdEnergy,"%g ",SysState.EnergyKineticFeedback);
#endif
/* comp */
for (i=0;i<6;i++)
{
fprintf(FdEnergy,"%g %g %g ",SysState.EnergyIntComp[i],SysState.EnergyPotComp[i], SysState.EnergyKinComp[i]);
#ifdef COOLING
fprintf(FdEnergy,"%g ",SysState.EnergyRadSphComp[i]);
#endif
#ifdef AGN_HEATING
fprintf(FdEnergy,"%g ",SysState.EnergyAGNHeatComp[i]);
#endif
#ifdef DISSIPATION_FORCES
fprintf(FdEnergy,"%g ",SysState.EnergyDissipationForcesComp[i]);
#endif
#ifdef INTEGRAL_CONSERVING_DISSIPATION
fprintf(FdEnergy,"%g ",SysState.EnergyICDissipationComp[i]);
#endif
#ifdef MULTIPHASE
fprintf(FdEnergy,"%g ",SysState.EnergyRadStickyComp[i]);
#endif
#ifdef FEEDBACK_WIND
fprintf(FdEnergy,"%g ",SysState.EnergyFeedbackWindComp[i]);
#endif
#ifdef BUBBLES
fprintf(FdEnergy,"%g ",SysState.EnergyBubblesComp[i]);
#endif
#ifdef CHIMIE_THERMAL_FEEDBACK
fprintf(FdEnergy,"%g ",SysState.EnergyThermalFeedbackComp[i]);
#endif
#ifdef CHIMIE_KINETIC_FEEDBACK
fprintf(FdEnergy,"%g ",SysState.EnergyKineticFeedbackComp[i]);
#endif
}
/* mass */
for (i=0;i<6;i++)
{
fprintf(FdEnergy,"%g ",SysState.MassComp[i]);
}
/* return */
fprintf(FdEnergy,"\n");
fflush(FdEnergy);
#ifdef SYSTEMSTATISTICS
/**************/
/* system */
/**************/
fprintf(FdSystem,"%g %g %g %g %g %g %g %g %g %g %g %g %g\n",
All.Time,
SysState.Momentum[0], SysState.Momentum[1], SysState.Momentum[2], SysState.Momentum[3],
SysState.AngMomentum[0], SysState.AngMomentum[1], SysState.AngMomentum[2], SysState.AngMomentum[3],
SysState.CenterOfMass[0], SysState.CenterOfMass[1], SysState.CenterOfMass[2], SysState.CenterOfMass[3]);
fflush(FdSystem);
#endif
}
#ifdef DETAILED_CPU
tend = second();
All.CPU_Residual += timediff(tstart, tend);
#endif
}
#endif
c4science · Help