statistics.c 12 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

/* Some standard headers. */
#include <string.h>

26
27
28
29
30
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif

31
32
33
34
35
36
37
/* This object's header. */
#include "statistics.h"

/* Local headers. */
#include "cooling.h"
#include "engine.h"
#include "error.h"
38
#include "gravity.h"
39
40
41
42
#include "hydro.h"
#include "threadpool.h"

/**
43
 * @brief Information required to compute the statistics in the mapper
44
45
 */
struct index_data {
46
  /*! The space we play with */
47
  const struct space *s;
48
49

  /*! The #statistics aggregator to fill */
50
51
52
  struct statistics *stats;
};

53
54
55
56
57
58
59
60
61
62
/**
 * @brief Adds the content of one #statistics aggregator to another one.
 *
 * Performs a += b;
 *
 * @param a The #statistics structure to update.
 * @param b The #statistics structure to add to a.
 */
void stats_add(struct statistics *a, const struct statistics *b) {

63
  /* Add everything */
64
65
  a->E_kin += b->E_kin;
  a->E_int += b->E_int;
66
67
  a->E_pot_self += b->E_pot_self;
  a->E_pot_ext += b->E_pot_ext;
68
69
70
71
72
73
74
75
76
  a->E_rad += b->E_rad;
  a->entropy += b->entropy;
  a->mass += b->mass;
  a->mom[0] += b->mom[0];
  a->mom[1] += b->mom[1];
  a->mom[2] += b->mom[2];
  a->ang_mom[0] += b->ang_mom[0];
  a->ang_mom[1] += b->ang_mom[1];
  a->ang_mom[2] += b->ang_mom[2];
77
78
79
  a->centre_of_mass[0] += b->centre_of_mass[0];
  a->centre_of_mass[1] += b->centre_of_mass[1];
  a->centre_of_mass[2] += b->centre_of_mass[2];
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
}

/**
 * @brief Initialises a statistics aggregator to a valid state.
 *
 * @param s The #statistics aggregator to initialise
 */
void stats_init(struct statistics *s) {

  /* Zero everything */
  bzero(s, sizeof(struct statistics));

  /* Set the lock */
  lock_init(&s->lock);
}

/**
97
 * @brief The #threadpool mapper function used to collect statistics for #part.
98
99
100
101
102
 *
 * @param map_data Pointer to the particles.
 * @param nr_parts The number of particles in this chunk
 * @param extra_data The #statistics aggregator.
 */
103
104
105
void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {

  /* Unpack the data */
106
  const struct index_data *data = (struct index_data *)extra_data;
107
  const struct space *s = data->s;
108
109
110
  const struct part *restrict parts = (struct part *)map_data;
  const struct xpart *restrict xparts =
      s->xparts + (ptrdiff_t)(parts - s->parts);
111
112
  const integertime_t ti_current = s->e->ti_current;
  const double timeBase = s->e->timeBase;
113
  const double time = s->e->time;
114
  struct statistics *const global_stats = data->stats;
115

116
117
118
119
  /* Required for external potential energy */
  const struct external_potential *potential = s->e->external_potential;
  const struct phys_const *phys_const = s->e->physical_constants;

120
121
  /* Local accumulator */
  struct statistics stats;
122
  stats_init(&stats);
123
124
125
126
127

  /* Loop over particles */
  for (int k = 0; k < nr_parts; k++) {

    /* Get the particle */
128
129
    const struct part *p = &parts[k];
    const struct xpart *xp = &xparts[k];
Matthieu Schaller's avatar
Matthieu Schaller committed
130
    const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL;
131

132
133
134
    /* Get useful time variables */
    const integertime_t ti_begin =
        get_integer_time_begin(ti_current, p->time_bin);
135
136
    const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
    const float dt = (ti_current - ((ti_begin + ti_end) / 2)) * timeBase;
137
138

    /* Get the total acceleration */
139
    float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
Matthieu Schaller's avatar
Matthieu Schaller committed
140
141
142
143
    if (gp != NULL) {
      a_tot[0] += gp->a_grav[0];
      a_tot[1] += gp->a_grav[1];
      a_tot[2] += gp->a_grav[2];
144
    }
145
146

    /* Extrapolate velocities to current time */
147
148
149
    const float v[3] = {xp->v_full[0] + a_tot[0] * dt,
                        xp->v_full[1] + a_tot[1] * dt,
                        xp->v_full[2] + a_tot[2] * dt};
150
151

    const float m = hydro_get_mass(p);
152
    const double x[3] = {p->x[0], p->x[1], p->x[2]};
153
154
    const float rho = hydro_get_density(p);
    const float u_int = hydro_get_internal_energy(p);
155

156
157
158
    /* Collect mass */
    stats.mass += m;

159
160
161
162
163
    /* Collect centre of mass */
    stats.centre_of_mass[0] += m * x[0];
    stats.centre_of_mass[1] += m * x[1];
    stats.centre_of_mass[2] += m * x[2];

164
165
166
167
168
169
170
171
172
173
174
175
    /* Collect momentum */
    stats.mom[0] += m * v[0];
    stats.mom[1] += m * v[1];
    stats.mom[2] += m * v[2];

    /* Collect angular momentum */
    stats.ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
    stats.ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
    stats.ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);

    /* Collect energies. */
    stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
176
    stats.E_int += m * u_int;
177
    stats.E_rad += cooling_get_radiated_energy(xp);
178
179
    if (gp != NULL) {
      stats.E_pot_self += m * gravity_get_potential(gp);
180
181
      stats.E_pot_ext += m * external_gravity_get_potential_energy(
                                 time, potential, phys_const, gp);
182
183
    }

184
    /* Collect entropy */
185
    stats.entropy += m * gas_entropy_from_internal_energy(rho, u_int);
186
187
188
  }

  /* Now write back to memory */
189
190
  if (lock_lock(&global_stats->lock) == 0) stats_add(global_stats, &stats);
  if (lock_unlock(&global_stats->lock) != 0) error("Failed to unlock stats.");
191
192
}

193
194
195
196
197
198
199
200
201
202
203
204
205
206
/**
 * @brief The #threadpool mapper function used to collect statistics for #gpart.
 *
 * @param map_data Pointer to the g-particles.
 * @param nr_gparts The number of g-particles in this chunk
 * @param extra_data The #statistics aggregator.
 */
void stats_collect_gpart_mapper(void *map_data, int nr_gparts,
                                void *extra_data) {

  /* Unpack the data */
  const struct index_data *data = (struct index_data *)extra_data;
  const struct space *s = data->s;
  const struct gpart *restrict gparts = (struct gpart *)map_data;
207
208
  const integertime_t ti_current = s->e->ti_current;
  const double timeBase = s->e->timeBase;
209
  const double time = s->e->time;
210
211
  struct statistics *const global_stats = data->stats;

212
213
214
215
  /* Required for external potential energy */
  const struct external_potential *potential = s->e->external_potential;
  const struct phys_const *phys_const = s->e->physical_constants;

216
217
218
219
220
221
222
223
224
225
226
227
228
229
  /* Local accumulator */
  struct statistics stats;
  stats_init(&stats);

  /* Loop over particles */
  for (int k = 0; k < nr_gparts; k++) {

    /* Get the particle */
    const struct gpart *gp = &gparts[k];

    /* If the g-particle has a counterpart, ignore it */
    if (gp->id_or_neg_offset < 0) continue;

    /* Get useful variables */
230
231
    const integertime_t ti_begin =
        get_integer_time_begin(ti_current, gp->time_bin);
232
233
    const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
    const float dt = (ti_current - ((ti_begin + ti_end) / 2)) * timeBase;
234
235

    /* Extrapolate velocities */
236
237
238
239
    const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt,
                        gp->v_full[1] + gp->a_grav[1] * dt,
                        gp->v_full[2] + gp->a_grav[2] * dt};

240
    const float m = gravity_get_mass(gp);
241
    const double x[3] = {gp->x[0], gp->x[1], gp->x[2]};
242
243
244
245

    /* Collect mass */
    stats.mass += m;

246
247
248
249
250
    /* Collect centre of mass */
    stats.centre_of_mass[0] += m * x[0];
    stats.centre_of_mass[1] += m * x[1];
    stats.centre_of_mass[2] += m * x[2];

251
252
253
254
255
256
257
258
259
260
261
262
    /* Collect momentum */
    stats.mom[0] += m * v[0];
    stats.mom[1] += m * v[1];
    stats.mom[2] += m * v[2];

    /* Collect angular momentum */
    stats.ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
    stats.ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
    stats.ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);

    /* Collect energies. */
    stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
263
    stats.E_pot_self += m * gravity_get_potential(gp);
264
265
    stats.E_pot_ext += m * external_gravity_get_potential_energy(
                               time, potential, phys_const, gp);
266
267
268
269
270
271
272
  }

  /* Now write back to memory */
  if (lock_lock(&global_stats->lock) == 0) stats_add(global_stats, &stats);
  if (lock_unlock(&global_stats->lock) != 0) error("Failed to unlock stats.");
}

273
274
275
276
277
278
/**
 * @brief Collect physical statistics over all particles in a #space.
 *
 * @param s The #space to collect from.
 * @param stats The #statistics aggregator to fill.
 */
279
280
281
282
283
void stats_collect(const struct space *s, struct statistics *stats) {

  /* Prepare the data */
  struct index_data extra_data;
  extra_data.s = s;
284
  extra_data.stats = stats;
285

286
287
288
  /* Run parallel collection of statistics for parts */
  if (s->nr_parts > 0)
    threadpool_map(&s->e->threadpool, stats_collect_part_mapper, s->parts,
289
                   s->nr_parts, sizeof(struct part), 0, &extra_data);
290
291
292
293

  /* Run parallel collection of statistics for gparts */
  if (s->nr_gparts > 0)
    threadpool_map(&s->e->threadpool, stats_collect_gpart_mapper, s->gparts,
294
                   s->nr_gparts, sizeof(struct gpart), 0, &extra_data);
295
}
296

297
298
299
300
301
302
303
304
305
306
307
308
/**
 * @brief Apply final opetations on the #stats.
 *
 * @param stats The #stats to work on.
 */
void stats_finalize(struct statistics *stats) {

  stats->centre_of_mass[0] /= stats->mass;
  stats->centre_of_mass[1] /= stats->mass;
  stats->centre_of_mass[2] /= stats->mass;
}

309
310
311
312
313
314
315
316
317
318
/**
 * @brief Prints the content of a #statistics aggregator to a file
 *
 * @param file File to write to.
 * @param stats The #statistics object to write to the file
 * @param time The current physical time.
 */
void stats_print_to_file(FILE *file, const struct statistics *stats,
                         double time) {

319
320
  const double E_pot = stats->E_pot_self + stats->E_pot_ext;
  const double E_tot = stats->E_kin + stats->E_int + E_pot;
321
322
323

  fprintf(file,
          " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e "
324
          "%14e %14e %14e %14e %14e %14e\n",
325
          time, stats->mass, E_tot, stats->E_kin, stats->E_int, E_pot,
Matthieu Schaller's avatar
Matthieu Schaller committed
326
327
          stats->E_pot_self, stats->E_pot_ext, stats->E_rad, stats->entropy,
          stats->mom[0], stats->mom[1], stats->mom[2], stats->ang_mom[0],
328
329
          stats->ang_mom[1], stats->ang_mom[2], stats->centre_of_mass[0],
          stats->centre_of_mass[1], stats->centre_of_mass[2]);
330
331
332
  fflush(file);
}

333
/* Extra stuff in MPI-land */
334
335
#ifdef WITH_MPI

336
337
338
/**
 * @brief MPI datatype corresponding to the #statistics structure.
 */
339
MPI_Datatype statistics_mpi_type;
340
341
342
343

/**
 * @brief MPI operator used for the reduction of #statistics structure.
 */
344
345
346
347
348
349
350
351
MPI_Op statistics_mpi_reduce_op;

/**
 * @brief MPI reduce operator for #statistics structures.
 */
void stats_add_MPI(void *in, void *inout, int *len, MPI_Datatype *datatype) {

  for (int i = 0; i < *len; ++i)
352
353
    stats_add(&((struct statistics *)inout)[0],
              &((const struct statistics *)in)[i]);
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
}

/**
 * @brief Registers MPI #statistics type and reduction function.
 */
void stats_create_MPI_type() {

  /* This is not the recommended way of doing this.
     One should define the structure field by field
     But as long as we don't do serialization via MPI-IO
     we don't really care.
     Also we would have to modify this function everytime something
     is added to the statistics structure. */
  if (MPI_Type_contiguous(sizeof(struct statistics) / sizeof(unsigned char),
                          MPI_BYTE, &statistics_mpi_type) != MPI_SUCCESS ||
      MPI_Type_commit(&statistics_mpi_type) != MPI_SUCCESS) {
    error("Failed to create MPI type for statistics.");
  }

  /* Create the reduction operation */
  MPI_Op_create(stats_add_MPI, 1, &statistics_mpi_reduce_op);
}
#endif