Collect the radiated energy of all particles when computing statistics and...

Collect the radiated energy of all particles when computing statistics and output the result in the log file.

src/cell.h

@@ -182,7 +182,7 @@ struct cell {
   double mom[3], ang_mom[3];
 
   /* Mass, potential, internal  and kinetic energy of particles in this cell. */
-  double mass, e_pot, e_int, e_kin, entropy;
+  double mass, e_pot, e_int, e_kin, e_rad, entropy;
 
   /* Number of particles updated in this cell. */
   int updated, g_updated;

src/cooling/const_du/cooling.h

@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   hydro_set_internal_energy(p, u_new);
 
   /* Store the radiated energy */
-  xp->cooling_data.radiated_energy += u_new - u_old;
+  xp->cooling_data.radiated_energy += hydro_get_mass(p) * (u_old - u_new);
 }
 
 /**
@@ -123,6 +123,20 @@ __attribute__((always_inline)) INLINE static void cooling_init_part(
   xp->cooling_data.radiated_energy = 0.f;
 }
 
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * In this simple example we jsut return the quantity accumulated in the
+ * #cooling_xpart_data of this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return xp->cooling_data.radiated_energy;
+}
+
 /**
  * @brief Initialises the cooling function properties from the parameter file
  *

src/cooling/const_du/cooling_struct.h

@@ -53,8 +53,7 @@ struct cooling_function_data {
  */
 struct cooling_xpart_data {
 
-  /*! Amount of energy radiated away by this particle since the start of the run
-   */
+  /*! Energy radiated away by this particle since the start of the run */
   float radiated_energy;
 };
 

src/cooling/const_lambda/cooling.h

@@ -116,7 +116,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   /*       u_old, du_dt, dt, du_dt * dt, u_new, u_new_test); */
 
   /* Store the radiated energy */
-  xp->cooling_data.radiated_energy += u_new - u_old;
+  xp->cooling_data.radiated_energy += hydro_get_mass(p) * (u_old - u_new);
 }
 
 /**
@@ -154,6 +154,17 @@ __attribute__((always_inline)) INLINE static void cooling_init_part(
   xp->cooling_data.radiated_energy = 0.f;
 }
 
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return xp->cooling_data.radiated_energy;
+}
+
 /**
  * @brief Initialises the cooling properties.
  *

src/cooling/const_lambda/cooling_struct.h

@@ -48,13 +48,12 @@ struct cooling_function_data {
   float cooling_tstep_mult;
 };
 
-
 /**
  * @brief Properties of the cooling stored in the particle data.
  */
 struct cooling_xpart_data {
 
-  /*! Amount of energy radiated away by this particle since the start of the run */
+  /*! Energy radiated away by this particle since the start of the run */
   float radiated_energy;
 };
 

src/cooling/none/cooling.h

@@ -83,6 +83,19 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
 __attribute__((always_inline)) INLINE static void cooling_init_part(
     const struct part* restrict p, struct xpart* restrict xp) {}
 
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * No cooling, so return 0.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+    const struct xpart* restrict xp) {
+
+  return 0.f;
+}
+
 /**
  * @brief Initialises the cooling properties.
  *

src/engine.c

@@ -2465,7 +2465,8 @@ void engine_print_stats(struct engine *e) {
   const ticks tic = getticks();
   const struct space *s = e->s;
 
-  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
+  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, e_rad = 0.0;
+  double entropy = 0.0, mass = 0.0;
   double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
 
   /* Collect the cell data. */
@@ -2476,6 +2477,7 @@ void engine_print_stats(struct engine *e) {
       e_kin += c->e_kin;
       e_int += c->e_int;
       e_pot += c->e_pot;
+      e_rad += c->e_rad;
       entropy += c->entropy;
       mom[0] += c->mom[0];
       mom[1] += c->mom[1];
@@ -2488,33 +2490,35 @@ void engine_print_stats(struct engine *e) {
 /* Aggregate the data from the different nodes. */
 #ifdef WITH_MPI
   {
-    double in[11] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
-    double out[11];
+    double in[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
+    double out[12];
     out[0] = e_kin;
     out[1] = e_int;
     out[2] = e_pot;
-    out[3] = mom[0];
-    out[4] = mom[1];
-    out[5] = mom[2];
-    out[6] = ang_mom[0];
-    out[7] = ang_mom[1];
-    out[8] = ang_mom[2];
-    out[9] = mass;
-    out[10] = entropy;
+    out[3] = e_rad;
+    out[4] = mom[0];
+    out[5] = mom[1];
+    out[6] = mom[2];
+    out[7] = ang_mom[0];
+    out[8] = ang_mom[1];
+    out[9] = ang_mom[2];
+    out[10] = mass;
+    out[11] = entropy;
     if (MPI_Reduce(out, in, 11, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) !=
         MPI_SUCCESS)
       error("Failed to aggregate stats.");
     e_kin = out[0];
     e_int = out[1];
     e_pot = out[2];
-    mom[0] = out[3];
-    mom[1] = out[4];
-    mom[2] = out[5];
-    ang_mom[0] = out[6];
-    ang_mom[1] = out[7];
-    ang_mom[2] = out[8];
-    mass = out[9];
-    entropy = out[10];
+    e_rad = out[3];
+    mom[0] = out[4];
+    mom[1] = out[5];
+    mom[2] = out[6];
+    ang_mom[0] = out[7];
+    ang_mom[1] = out[8];
+    ang_mom[2] = out[9];
+    mass = out[10];
+    entropy = out[11];
   }
 #endif
 
@@ -2522,11 +2526,11 @@ void engine_print_stats(struct engine *e) {
 
   /* Print info */
   if (e->nodeID == 0) {
-    fprintf(
-        e->file_stats,
-        " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
-        e->time, mass, e_tot, e_kin, e_int, e_pot, entropy, mom[0], mom[1],
-        mom[2], ang_mom[0], ang_mom[1], ang_mom[2]);
+    fprintf(e->file_stats,
+            " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e "
+            "%14e\n",
+            e->time, mass, e_tot, e_kin, e_int, e_pot, e_rad, entropy, mom[0],
+            mom[1], mom[2], ang_mom[0], ang_mom[1], ang_mom[2]);
     fflush(e->file_stats);
   }
 
@@ -3358,11 +3362,11 @@ void engine_init(struct engine *e, struct space *s,
                                 engine_default_energy_file_name);
     sprintf(energyfileName + strlen(energyfileName), ".txt");
     e->file_stats = fopen(energyfileName, "w");
-    fprintf(
-        e->file_stats,
-        "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
-        "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "Entropy", "p_x",
-        "p_y", "p_z", "ang_x", "ang_y", "ang_z");
+    fprintf(e->file_stats,
+            "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
+            "%14s\n",
+            "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "E_rad",
+            "Entropy", "p_x", "p_y", "p_z", "ang_x", "ang_y", "ang_z");
     fflush(e->file_stats);
 
     char timestepsfileName[200] = "";

src/runner.c

@@ -734,7 +734,8 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
   const double dt = (ti_current - ti_old) * timeBase;
 
   float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
-  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
+  double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, e_rad = 0.0;
+  double entropy = 0.0, mass = 0.0;
   double mom[3] = {0.0, 0.0, 0.0};
   double ang_mom[3] = {0.0, 0.0, 0.0};
 
@@ -806,6 +807,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
       e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
       e_pot += 0.;
       e_int += m * hydro_get_internal_energy(p, half_dt);
+      e_rad += cooling_get_radiated_energy(xp);
 
       /* Collect entropy */
       entropy += m * hydro_get_entropy(p, half_dt);
@@ -832,6 +834,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
         e_kin += cp->e_kin;
         e_int += cp->e_int;
         e_pot += cp->e_pot;
+        e_rad += cp->e_rad;
         entropy += cp->entropy;
         mom[0] += cp->mom[0];
         mom[1] += cp->mom[1];
@@ -849,6 +852,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
   c->e_kin = e_kin;
   c->e_int = e_int;
   c->e_pot = e_pot;
+  c->e_rad = e_rad;
   c->entropy = entropy;
   c->mom[0] = mom[0];
   c->mom[1] = mom[1];