From f0e38086f37b8ee635a196421c0ea966a39d0cbb Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <vandenbroucke@strw.leidenuniv.nl>
Date: Wed, 31 Aug 2022 13:27:10 +0000
Subject: [PATCH] Ghost iteration statistics
---
configure.ac | 15 +
doc/RTD/source/AnalysisTools/index.rst | 47 +++
src/Makefile.am | 2 +
src/cell.h | 3 +
src/engine.c | 5 +
src/ghost_stats.c | 111 ++++++
src/ghost_stats.h | 500 +++++++++++++++++++++++++
src/runner_ghost.c | 25 ++
swift.c | 2 +
tools/plot_ghost_stats.py | 229 +++++++++++
10 files changed, 939 insertions(+)
create mode 100644 src/ghost_stats.c
create mode 100644 src/ghost_stats.h
create mode 100644 tools/plot_ghost_stats.py
diff --git a/configure.ac b/configure.ac
index 04c82357c0..9357d35df0 100644
--- a/configure.ac
+++ b/configure.ac
@@ -405,6 +405,20 @@ elif test "$stars_density_checks" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_STARS_DENSITY_CHECKS], [$enableval] ,[Enable stars density brute-force checks])
fi
+# Check if ghost statistics are enabled
+AC_ARG_ENABLE([ghost-statistics],
+ [AS_HELP_STRING([--enable-ghost-statistics],
+ [Gather statistics about the ghost iterations for hydro, stars and black holes in N bins @<:@N@:>@]
+ )],
+ [ghost_stats="$enableval"],
+ [ghost_stats="no"]
+)
+if test "$ghost_stats" == "yes"; then
+ AC_MSG_ERROR(Need to specify the number of bins when using --enable-ghost-statistics!)
+elif test "$ghost_stats" != "no"; then
+ AC_DEFINE_UNQUOTED([SWIFT_GHOST_STATS], [$enableval] ,[Enable ghost statistics for hydro, stars and black holes])
+fi
+
# Check whether we want to switch on glass making
AC_ARG_ENABLE([glass-making],
[AS_HELP_STRING([--enable-glass-making],
@@ -2984,6 +2998,7 @@ AC_MSG_RESULT([
Boundary particles : $boundary_particles
Fixed boundary particles : $fixed_boundary_particles
Planetary fixed entropy : $planetary_fixed_entropy
+ Ghost statistics : $ghost_stats
Continuous Sim. Data Stream : $with_csds
diff --git a/doc/RTD/source/AnalysisTools/index.rst b/doc/RTD/source/AnalysisTools/index.rst
index f86a793072..15922608cf 100644
--- a/doc/RTD/source/AnalysisTools/index.rst
+++ b/doc/RTD/source/AnalysisTools/index.rst
@@ -235,3 +235,50 @@ The ``.dump<.rank>`` files once seen are deleted, so dumping can be done more
than once. For a non-MPI run the file is simply called ``.dump``, note for MPI
you need to create one file per rank, so ``.dump.0``, ``.dump.1`` and so on.
+
+Neighbour search statistics
+---------------------------
+
+One of the core algorithms in SWIFT is an iterative neighbour search
+whereby we try to find an appropriate radius around a particle's
+position so that the weighted sum over neighbouring particles within
+that radius is equal to some target value. The most obvious example of
+this iterative neighbour search is the SPH density loop, but various
+sub-grid models employ a very similar iterative neighbour search. The
+computational cost of this iterative search is significantly affected by
+the number of iterations that is required, and it can therefore be
+useful to analyse the progression of the iterative scheme in detail.
+
+When configured with ``--enable-ghost-statistics=X``, SWIFT will be
+compiled with additional diagnostics that statistically track the number
+of iterations required to find a converged answer. Here, ``X`` is a
+fixed number of bins to use to collect the required statistics
+(``ghost`` refers to the fact that the iterations take place inside the
+ghost tasks). In practice, this means that every cell in the SWIFT tree
+will be equipped with an additional ``struct`` containing three sets of
+``X`` bins (one set for each iterative neighbour loop: hydro, stellar
+feedback, AGN feedback). For each bin ``i``, we store the number of
+particles that required updating during iteration ``i``, the number of
+particles that could not find a single neighbouring particle, the
+minimum and maximum smoothing length of all particles that required
+updating, and the sum of all their search radii and all their search
+radii squared. This allows us to calculate the upper and lower limits,
+as well as the mean and standard deviation on the search radius for each
+iteration and for each cell. Note that there could be more iterations
+required than the number of bins ``X``; in this case the additional
+iterations will be accumulated in the final bin. At the end of each time
+step, a text file is produced (one per MPI rank) that contains the
+information for all cells that had any relevant activity. This text file
+is named ``ghost_stats_ssss_rrrr.txt``, where ``ssss`` is the step
+counter for that time step and ``rrrr`` is the MPI rank.
+
+The script ``tools/plot_ghost_stats.py`` takes one or multiple
+``ghost_stats.txt`` files and computes global statistics for all the
+cells in those files. The script also takes the name of an output file
+where it will save those statistics as a set of plots, and an optional
+label that will be displayed as the title of the plots. Note that there
+are no restrictions on the number of input files or how they relate;
+different files could represent different MPI ranks, but also different
+time steps or even different simulations (which would make little
+sense). It is up to the user to make sure that the input is actually
+relevant.
diff --git a/src/Makefile.am b/src/Makefile.am
index 3177bff00e..3dbb2cfbf4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -73,6 +73,7 @@ include_HEADERS += lightcone/lightcone_crossing.h lightcone/lightcone_array.h li
include_HEADERS += lightcone/lightcone_map_types.h lightcone/projected_kernel.h lightcone/lightcone_shell.h
include_HEADERS += lightcone/healpix_util.h lightcone/pixel_index.h
include_HEADERS += power_spectrum.h
+include_HEADERS += ghost_stats.h
# source files for EAGLE extra I/O
EAGLE_EXTRA_IO_SOURCES=
@@ -177,6 +178,7 @@ AM_SOURCES += lightcone/lightcone.c lightcone/lightcone_particle_io.c lightcone/
AM_SOURCES += lightcone/healpix_util.c lightcone/lightcone_array.c lightcone/lightcone_map.c
AM_SOURCES += lightcone/lightcone_map_types.c lightcone/projected_kernel.c lightcone/lightcone_shell.c
AM_SOURCES += power_spectrum.c
+AM_SOURCES += ghost_stats.c
AM_SOURCES += $(EAGLE_EXTRA_IO_SOURCES)
AM_SOURCES += $(QLA_COOLING_SOURCES) $(QLA_EAGLE_COOLING_SOURCES)
AM_SOURCES += $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES)
diff --git a/src/cell.h b/src/cell.h
index 6fd8b7d95b..da6c610c73 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -38,6 +38,7 @@
#include "cell_hydro.h"
#include "cell_sinks.h"
#include "cell_stars.h"
+#include "ghost_stats.h"
#include "kernel_hydro.h"
#include "multipole_struct.h"
#include "part.h"
@@ -463,6 +464,8 @@ struct cell {
char subtasks_executed[task_type_count];
#endif
+ struct ghost_stats ghost_statistics;
+
} SWIFT_STRUCT_ALIGN;
/* Convert cell location to ID. */
diff --git a/src/engine.c b/src/engine.c
index 86faa07863..bbae68c777 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2398,6 +2398,11 @@ int engine_step(struct engine *e) {
e->step % e->sched.frequency_task_levels == 0)
scheduler_write_task_level(&e->sched, e->step);
+ /* we have to reset the ghost histograms here and not in engine_launch,
+ because engine_launch is re-used for the limiter and sync (and we don't
+ want to lose the data from the tasks) */
+ space_reset_ghost_histograms(e->s);
+
/* Start all the tasks. */
TIMER_TIC;
engine_launch(e, "tasks");
diff --git a/src/ghost_stats.c b/src/ghost_stats.c
new file mode 100644
index 0000000000..883579d9f5
--- /dev/null
+++ b/src/ghost_stats.c
@@ -0,0 +1,111 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#include "ghost_stats.h"
+
+#include "cell.h"
+#include "space.h"
+
+#ifdef SWIFT_GHOST_STATS
+
+/**
+ * @brief Resets all the individual cell ghost counters to 0.
+ *
+ * @param c The #cell to reset.
+ */
+void cell_reset_ghost_histograms(struct cell *c) {
+
+ ghost_stats_reset_entries(&c->ghost_statistics);
+
+ /* recurse */
+ for (int k = 0; k < 8; ++k)
+ if (c->progeny[k] != NULL) cell_reset_ghost_histograms(c->progeny[k]);
+}
+
+/**
+ * @brief Write the ghost statistics for the given cell to the given file,
+ * using the given cell ID.
+ *
+ * This function recurses into the progenies. They all get different cell IDs
+ * by shifting the parent cell ID by 3 bits. Note that this does not
+ * guarantee unique cell IDs (but it is good enough to distinguish parents
+ * from their progeny).
+ *
+ * @param f File to write to.
+ * @param c Cell to write.
+ * @param cellID ID used to distinguish this cell from its progeny and
+ * immediate neigbours (not guaranteed to be unique!).
+ */
+void cell_write_ghost_stats(FILE *f, const struct cell *c,
+ const long long cellID) {
+ if (c == NULL) return;
+
+ ghost_stats_write_cell_stats(f, &c->ghost_statistics, cellID);
+
+ /* Write children */
+ const long long pID = cellID << 3;
+ for (int i = 0; i < 8; i++) {
+ cell_write_ghost_stats(f, c->progeny[i], pID + i);
+ }
+}
+
+/**
+ * @brief Reset the ghost histograms for all top level cells in the space.
+ *
+ * @param s Space.
+ */
+void space_reset_ghost_histograms(struct space *s) {
+ for (int i = 0; i < s->nr_cells; ++i) {
+ cell_reset_ghost_histograms(&s->cells_top[i]);
+ }
+}
+
+/**
+ * @brief Write the ghost statistics for all the top level cells in the space
+ * to a file.
+ *
+ * The file name is composed as follows:
+ * ghost_stats_%04i_%04i.txt
+ * where the first counter is an argument to this function (intended to be the
+ * step counter), and the second counter is the rank that does the write.
+ *
+ * @param s Space.
+ * @param j First counter in the output file name.
+ */
+void space_write_ghost_stats(const struct space *s, int j) {
+ /* Open file */
+ char filename[200];
+ sprintf(filename, "ghost_stats_%04i_%04i.txt", j, engine_rank);
+ FILE *f = fopen(filename, "w");
+ if (f == NULL) error("Error opening ghost stats file.");
+
+ /* Write header */
+ ghost_stats_write_header(f);
+
+ /* Write all the top level cells (and their children) */
+ for (int i = 0; i < s->nr_cells; i++) {
+ struct cell *c = &s->cells_top[i];
+ if (c->nodeID == engine_rank) cell_write_ghost_stats(f, c, i);
+ }
+
+ /* Cleanup */
+ fclose(f);
+}
+
+#endif
diff --git a/src/ghost_stats.h b/src/ghost_stats.h
new file mode 100644
index 0000000000..23fbcc288c
--- /dev/null
+++ b/src/ghost_stats.h
@@ -0,0 +1,500 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2021 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+#ifndef SWIFT_GHOST_STATS_H
+#define SWIFT_GHOST_STATS_H
+
+/* Config parameters. */
+#include "minmax.h"
+#include "part.h"
+
+#include <config.h>
+#include <float.h>
+#include <stdio.h>
+
+#ifdef SWIFT_GHOST_STATS
+
+/*****************************
+ * ghost_stats functionality *
+ *****************************/
+
+/**
+ * @brief Entry for a single iteration in the ghost statistics table.
+ */
+struct ghost_stats_entry {
+ /*! Number of particles processed during the iteration. */
+ int count;
+ /*! Number of particles in the iteration that had no neighbours. */
+ int count_no_ngb;
+ /*! Minimum initial smoothing length of the particles at the start of the
+ * iteration. */
+ float hmin;
+ /*! Maximum initial smoothing length of the particles at the start of the
+ * iteration. */
+ float hmax;
+ /*! Sum of the initial smoothing lengths, useful to compute averages in
+ * post-processing. */
+ double hsum;
+ /*! Sum of the initial smoothing lengths squared, useful to compute variances
+ * and standard deviations in post-processing. */
+ double hsum2;
+};
+
+/**
+ * @brief Ghost statistics stored in a cell.
+ */
+struct ghost_stats {
+ /* Hydro ghost statistics. */
+ struct ghost_stats_entry hydro[SWIFT_GHOST_STATS + 1];
+ /* Stars ghost statistics. */
+ struct ghost_stats_entry stars[SWIFT_GHOST_STATS + 1];
+ /* Black holes ghost statistics. */
+ struct ghost_stats_entry black_holes[SWIFT_GHOST_STATS + 1];
+};
+
+/* ghost_stats_entry struct functions */
+
+/**
+ * @brief Reset the given ghost stats bin.
+ *
+ * @param bin Ghost stats bin to reset.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_reset_entry(
+ struct ghost_stats_entry *restrict bin) {
+
+ bin->count = 0;
+ bin->count_no_ngb = 0;
+ bin->hmin = FLT_MAX;
+ bin->hmax = 0.0f;
+ bin->hsum = 0.;
+ bin->hsum2 = 0.;
+}
+
+/**
+ * @brief Write the given ghost stats bin to the given file.
+ *
+ * @param f File to write to.
+ * @param bin Ghost stats bin to write.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_print_entry(
+ FILE *f, const struct ghost_stats_entry *restrict bin) {
+
+ if (bin->count > 0) {
+ fprintf(f, "\t%i\t%i\t%g\t%g\t%g\t%g", bin->count, bin->count_no_ngb,
+ bin->hmin, bin->hmax, bin->hsum, bin->hsum2);
+ } else {
+ fprintf(f, "\t0\t0\t0\t0\t0\t0");
+ }
+}
+
+/* ghost_stats struct functions */
+
+/**
+ * @brief Reset all the entries in the ghost_stats struct.
+ *
+ * @param gstats Ghost stats struct.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_reset_entries(
+ struct ghost_stats *restrict gstats) {
+
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_reset_entry(&gstats->hydro[b]);
+ }
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_reset_entry(&gstats->stars[b]);
+ }
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_reset_entry(&gstats->black_holes[b]);
+ }
+}
+
+/**
+ * @brief Account for the star particles that are still under consideration at
+ * the start of a ghost decision loop (so after the neighbour loop but before
+ * the smoothing length is updated).
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Iteration number in the high-level ghost iteration
+ * scheme.
+ * @param scount Number of star particles still under consideration (smoothing
+ * length has not converged yet or will do so during this iteration).
+ * @param sparts Star particle array.
+ * @param sid Indices of active star particles in the array.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_account_for_stars(
+ struct ghost_stats *restrict gstats, int iteration_number, int scount,
+ struct spart *restrict sparts, int *sid) {
+
+ /* accumulate in highest bin, so we can spot out-of-bounds values */
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict sbin = &gstats->stars[binidx];
+ sbin->count += scount;
+ for (int i = 0; i < scount; i++) {
+ const float hi = sparts[sid[i]].h;
+ sbin->hmin = min(sbin->hmin, hi);
+ sbin->hmax = max(sbin->hmax, hi);
+ sbin->hsum += hi;
+ sbin->hsum2 += hi * hi;
+ }
+}
+
+/**
+ * @brief Account for the properties of the a converged star particle.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param sp Star particle that has a converged smoothing length.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_converged_star(
+ struct ghost_stats *restrict gstats, struct spart *restrict sp) {
+
+ struct ghost_stats_entry *restrict sbin = &gstats->stars[SWIFT_GHOST_STATS];
+ ++sbin->count;
+ const float hi = sp->h;
+ sbin->hmin = min(sbin->hmin, hi);
+ sbin->hmax = max(sbin->hmax, hi);
+ sbin->hsum += hi;
+ sbin->hsum2 += hi * hi;
+}
+
+/**
+ * @brief Register the occurrence of a "no neighbour" event during the current
+ * star iteration step.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Number of the current iteration.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_star_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {
+
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict sbin = &gstats->stars[binidx];
+ ++sbin->count_no_ngb;
+}
+
+/**
+ * @brief Register the occurrence of a converged star particle without
+ * neighbours.
+ *
+ * @param gstats Ghost stats struct to update.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_star_converged(struct ghost_stats *restrict gstats) {
+
+ struct ghost_stats_entry *restrict sbin = &gstats->stars[SWIFT_GHOST_STATS];
+ ++sbin->count_no_ngb;
+}
+
+/**
+ * @brief Account for the black hole particles that are still under
+ * consideration at the start of a ghost decision loop (so after the neighbour
+ * loop but before the smoothing length is updated).
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Iteration number in the high-level ghost iteration
+ * scheme.
+ * @param bcount Number of black hole particles still under consideration
+ * (smoothing length has not converged yet or will do so during this iteration).
+ * @param bparts Black hole particle array.
+ * @param sid Indices of active black hole particles in the array.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_account_for_black_holes(struct ghost_stats *restrict gstats,
+ int iteration_number, int bcount,
+ struct bpart *restrict bparts, int *sid) {
+
+ /* accumulate in highest bin, so we can spot out-of-bounds values */
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict bbin = &gstats->black_holes[binidx];
+ bbin->count += bcount;
+ for (int i = 0; i < bcount; i++) {
+ const float hi = bparts[sid[i]].h;
+ bbin->hmin = min(bbin->hmin, hi);
+ bbin->hmax = max(bbin->hmax, hi);
+ bbin->hsum += hi;
+ bbin->hsum2 += hi * hi;
+ }
+}
+
+/**
+ * @brief Account for the properties of the a converged black hole particle.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param bp Black hole particle that has a converged smoothing length.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_converged_black_hole(struct ghost_stats *restrict gstats,
+ struct bpart *restrict bp) {
+
+ struct ghost_stats_entry *restrict bbin =
+ &gstats->black_holes[SWIFT_GHOST_STATS];
+ ++bbin->count;
+ const float hi = bp->h;
+ bbin->hmin = min(bbin->hmin, hi);
+ bbin->hmax = max(bbin->hmax, hi);
+ bbin->hsum += hi;
+ bbin->hsum2 += hi * hi;
+}
+
+/**
+ * @brief Register the occurrence of a "no neighbour" event during the current
+ * black hole iteration step.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Number of the current iteration.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_black_hole_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {
+
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict bbin = &gstats->black_holes[binidx];
+ ++bbin->count_no_ngb;
+}
+
+/**
+ * @brief Register the occurrence of a converged black hole particle without
+ * neighbours.
+ *
+ * @param gstats Ghost stats struct to update.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_black_hole_converged(struct ghost_stats *restrict gstats) {
+
+ struct ghost_stats_entry *restrict bbin =
+ &gstats->black_holes[SWIFT_GHOST_STATS];
+ ++bbin->count_no_ngb;
+}
+
+/**
+ * @brief Account for the gas particles that are still under consideration at
+ * the start of a ghost decision loop (so after the neighbour loop but before
+ * the smoothing length is updated).
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Iteration number in the high-level ghost iteration
+ * scheme.
+ * @param count Number of gas particles still under consideration (smoothing
+ * length has not converged yet or will do so during this iteration).
+ * @param parts Gas particle array.
+ * @param sid Indices of active gas particles in the array.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_account_for_hydro(
+ struct ghost_stats *restrict gstats, int iteration_number, int count,
+ struct part *restrict parts, int *pid) {
+
+ /* accumulate in highest bin, so we can spot out-of-bounds values */
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict hbin = &gstats->hydro[binidx];
+ hbin->count += count;
+ for (int i = 0; i < count; i++) {
+ const float hi = parts[pid[i]].h;
+ hbin->hmin = min(hbin->hmin, hi);
+ hbin->hmax = max(hbin->hmax, hi);
+ hbin->hsum += hi;
+ hbin->hsum2 += hi * hi;
+ }
+}
+
+/**
+ * @brief Account for the properties of the a converged gas particle.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param p Gas particle that has a converged smoothing length.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_converged_hydro(
+ struct ghost_stats *restrict gstats, struct part *restrict p) {
+
+ struct ghost_stats_entry *restrict hbin = &gstats->hydro[SWIFT_GHOST_STATS];
+ ++hbin->count;
+ const float hi = p->h;
+ hbin->hmin = min(hbin->hmin, hi);
+ hbin->hmax = max(hbin->hmax, hi);
+ hbin->hsum += hi;
+ hbin->hsum2 += hi * hi;
+}
+
+/**
+ * @brief Register the occurrence of a "no neighbour" event during the current
+ * hydro iteration step.
+ *
+ * @param gstats Ghost stats struct to update.
+ * @param iteration_number Number of the current iteration.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_hydro_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {
+
+ int binidx = min(iteration_number, SWIFT_GHOST_STATS - 1);
+ struct ghost_stats_entry *restrict hbin = &gstats->hydro[binidx];
+ ++hbin->count_no_ngb;
+}
+
+/**
+ * @brief Register the occurrence of a converged gas particle without
+ * neighbours.
+ *
+ * @param gstats Ghost stats struct to update.
+ */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_hydro_converged(struct ghost_stats *restrict gstats) {
+
+ struct ghost_stats_entry *restrict hbin = &gstats->hydro[SWIFT_GHOST_STATS];
+ ++hbin->count_no_ngb;
+}
+
+/**
+ * @brief Write the header of a ghost statistics file.
+ *
+ * @param f File to write to.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_write_header(
+ FILE *f) {
+
+ fprintf(f, "# Ghost statistics\n");
+ fprintf(f, "# Values listed in blocks per particle type\n");
+ fprintf(f, "# Order of types: hydro, stars, black holes\n");
+ fprintf(f, "# Number of blocks per type: %i\n", SWIFT_GHOST_STATS + 1);
+ fprintf(f, "# Number of values per block: 6\n");
+ fprintf(f, "# Last block contains converged values\n");
+ fprintf(f, "# Fields per block:\n");
+ fprintf(f, "# - count: i4\n");
+ fprintf(f, "# - no neighbour count: i4\n");
+ fprintf(f, "# - min h: f4\n");
+ fprintf(f, "# - max h: f4\n");
+ fprintf(f, "# - sum h: f8\n");
+ fprintf(f, "# - sum h^2: f8\n");
+ fprintf(f, "# First column is cellID\n");
+ fprintf(f, "# Cells with no values are omitted\n");
+}
+
+/**
+ * @brief Write the ghost statistics for the given cell to the given file.
+ *
+ * @param f File to write to.
+ * @param gstats Ghost statistics to write.
+ * @param cellID Cell ID to use to identify the cell.
+ */
+__attribute__((always_inline)) INLINE static void ghost_stats_write_cell_stats(
+ FILE *f, const struct ghost_stats *restrict gstats,
+ const long long cellID) {
+
+ if (gstats->hydro[0].count + gstats->stars[0].count > 0) {
+ fprintf(f, "%lld", cellID);
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_print_entry(f, &gstats->hydro[b]);
+ }
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_print_entry(f, &gstats->stars[b]);
+ }
+ for (int b = 0; b < SWIFT_GHOST_STATS + 1; ++b) {
+ ghost_stats_print_entry(f, &gstats->black_holes[b]);
+ }
+ fprintf(f, "\n");
+ }
+}
+
+/******************
+ * cell interface *
+ ******************/
+
+struct cell;
+
+void cell_reset_ghost_histograms(struct cell *c);
+void cell_write_ghost_stats(FILE *f, const struct cell *c,
+ const long long cellID);
+/*******************
+ * space interface *
+ *******************/
+
+struct space;
+
+void space_reset_ghost_histograms(struct space *s);
+void space_write_ghost_stats(const struct space *s, int j);
+
+#else
+
+struct ghost_stats {};
+
+/* stars */
+__attribute__((always_inline)) INLINE static void ghost_stats_account_for_stars(
+ struct ghost_stats *restrict gstats, int iteration_number, int scount,
+ struct spart *restrict sparts, int *sid) {}
+
+__attribute__((always_inline)) INLINE static void ghost_stats_converged_star(
+ struct ghost_stats *restrict gstats, struct spart *restrict sp) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_star_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_star_converged(struct ghost_stats *restrict gstats) {}
+
+/* black holes */
+__attribute__((always_inline)) INLINE static void
+ghost_stats_account_for_black_holes(struct ghost_stats *restrict gstats,
+ int iteration_number, int bcount,
+ struct bpart *restrict bparts, int *sid) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_converged_black_hole(struct ghost_stats *restrict gstats,
+ struct bpart *restrict bp) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_black_hole_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_black_hole_converged(struct ghost_stats *restrict gstats) {}
+
+/* hydro */
+__attribute__((always_inline)) INLINE static void ghost_stats_account_for_hydro(
+ struct ghost_stats *restrict gstats, int iteration_number, int count,
+ struct part *restrict parts, int *pid) {}
+
+__attribute__((always_inline)) INLINE static void ghost_stats_converged_hydro(
+ struct ghost_stats *restrict gstats, struct part *restrict p) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_hydro_iteration(struct ghost_stats *restrict gstats,
+ int iteration_number) {}
+
+__attribute__((always_inline)) INLINE static void
+ghost_stats_no_ngb_hydro_converged(struct ghost_stats *restrict gstats) {}
+
+/// cell interface
+
+struct cell;
+
+__attribute__((always_inline)) INLINE static void cell_reset_ghost_histograms(
+ struct cell *c) {}
+
+/// space interface
+
+struct space;
+
+__attribute__((always_inline)) INLINE static void space_reset_ghost_histograms(
+ struct space *s) {}
+
+__attribute__((always_inline)) INLINE static void space_write_ghost_stats(
+ const struct space *s, int j) {}
+#endif
+
+#endif /* SWIFT_GHOST_STATS */
diff --git a/src/runner_ghost.c b/src/runner_ghost.c
index 034c4c5af5..8750bd5089 100644
--- a/src/runner_ghost.c
+++ b/src/runner_ghost.c
@@ -145,6 +145,9 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
for (int num_reruns = 0; scount > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
+ ghost_stats_account_for_stars(&c->ghost_statistics, num_reruns, scount,
+ sparts, sid);
+
/* Reset the redo-count. */
redo = 0;
@@ -173,6 +176,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
if (sp->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */
+ ghost_stats_no_ngb_star_iteration(&c->ghost_statistics, num_reruns);
+
/* Flag that there were no neighbours */
has_no_neighbours = 1;
@@ -364,6 +369,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Do some damage control if no neighbours at all were found */
if (has_no_neighbours) {
+ ghost_stats_no_ngb_star_converged(&c->ghost_statistics);
stars_spart_has_no_neighbours(sp, cosmo);
rt_spart_has_no_neighbours(sp);
}
@@ -381,6 +387,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
h_max = max(h_max, sp->h);
h_max_active = max(h_max_active, sp->h);
+ ghost_stats_converged_star(&c->ghost_statistics, sp);
+
stars_reset_feedback(sp);
/* Only do feedback if stars have a reasonable birth time */
@@ -632,6 +640,9 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
for (int num_reruns = 0; bcount > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
+ ghost_stats_account_for_black_holes(&c->ghost_statistics, num_reruns,
+ bcount, bparts, sid);
+
/* Reset the redo-count. */
redo = 0;
@@ -658,6 +669,9 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
if (bp->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */
+ ghost_stats_no_ngb_black_hole_iteration(&c->ghost_statistics,
+ num_reruns);
+
/* Flag that there were no neighbours */
has_no_neighbours = 1;
@@ -774,6 +788,7 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
/* Do some damage control if no neighbours at all were found */
if (has_no_neighbours) {
+ ghost_stats_no_ngb_black_hole_converged(&c->ghost_statistics);
black_holes_bpart_has_no_neighbours(bp, cosmo);
}
@@ -786,6 +801,8 @@ void runner_do_black_holes_density_ghost(struct runner *r, struct cell *c,
/* We now have a particle whose smoothing length has converged */
+ ghost_stats_converged_black_hole(&c->ghost_statistics, bp);
+
black_holes_reset_feedback(bp);
/* Check if h_max has increased */
@@ -1135,6 +1152,9 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
for (int num_reruns = 0; count > 0 && num_reruns < max_smoothing_iter;
num_reruns++) {
+ ghost_stats_account_for_hydro(&c->ghost_statistics, num_reruns, count,
+ parts, pid);
+
/* Reset the redo-count. */
redo = 0;
@@ -1161,6 +1181,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
if (p->density.wcount < 1.e-5 * kernel_root) { /* No neighbours case */
+ ghost_stats_no_ngb_hydro_iteration(&c->ghost_statistics, num_reruns);
+
/* Flag that there were no neighbours */
has_no_neighbours = 1;
@@ -1366,6 +1388,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Do some damage control if no neighbours at all were found */
if (has_no_neighbours) {
+ ghost_stats_no_ngb_hydro_converged(&c->ghost_statistics);
hydro_part_has_no_neighbours(p, xp, cosmo);
mhd_part_has_no_neighbours(p, xp, cosmo);
chemistry_part_has_no_neighbours(p, xp, chemistry, cosmo);
@@ -1388,6 +1411,8 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
h_max = max(h_max, p->h);
h_max_active = max(h_max_active, p->h);
+ ghost_stats_converged_hydro(&c->ghost_statistics, p);
+
#ifdef EXTRA_HYDRO_LOOP
/* As of here, particle gradient variables will be set. */
diff --git a/swift.c b/swift.c
index 7be7fb95a3..94740abbe5 100644
--- a/swift.c
+++ b/swift.c
@@ -1700,6 +1700,8 @@ int main(int argc, char *argv[]) {
}
#endif
+ space_write_ghost_stats(e.s, j + 1);
+
/* Dump memory use report if collected. */
#ifdef SWIFT_MEMUSE_REPORTS
{
diff --git a/tools/plot_ghost_stats.py b/tools/plot_ghost_stats.py
new file mode 100644
index 0000000000..6e88043076
--- /dev/null
+++ b/tools/plot_ghost_stats.py
@@ -0,0 +1,229 @@
+#!/usr/bin/env python
+#
+# Usage:
+# python plot_ghost_stats.py ghost_stats_file_1 [ghost_stats_file_2...]
+# output_image.png
+#
+# Description:
+# Plot the ghost statistics contained in the given input files. For each
+# particle type (hydro, stars, black holes), three plots are generated:
+# 1. A histogram of the number of active particles as a function of iteration
+# number.
+# 2. The distribution of h values for active particles as a function of the
+# iteration number, compared with the final (converged) values.
+# 3. The number of particles that are treated as having no neighbours as a
+# function of the iteration number.
+
+import numpy as np
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as pl
+import argparse
+
+argparser = argparse.ArgumentParser()
+argparser.add_argument("input", nargs="+", action="store")
+argparser.add_argument("output", action="store")
+argparser.add_argument("--label", "-l", action="store", default=None)
+args = argparser.parse_args()
+
+data = None
+nbin = None
+nval = None
+for input in args.input:
+ # read the header block to figure out how many blocks per particle type
+ # we have
+ with open(input, "r") as ifile:
+ header = ifile.readlines()[:15]
+ this_nbin = int(header[3].split()[-1])
+ this_nval = int(header[4].split()[-1])
+ if nbin is None:
+ nbin = this_nbin
+ else:
+ if nbin != this_nbin:
+ print("Number of bins in files does not match!")
+ exit(1)
+ if nval is None:
+ nval = this_nval
+ else:
+ if nval != this_nval:
+ print("Number of values per bin in files does not match!")
+ exit(1)
+
+ tdata = np.loadtxt(input, ndmin=2)
+ if len(tdata) > 0:
+ if data is None:
+ data = np.array(tdata)
+ else:
+ data = np.append(data, tdata, axis=0)
+
+# omit meaningless first column
+data = data[:, 1:]
+
+hdata = data[:, : nval * nbin]
+sdata = data[:, nval * nbin : 2 * nval * nbin]
+bdata = data[:, 2 * nval * nbin :]
+
+
+def get_total_stats(pdata):
+ counts = pdata[:, ::nval]
+ ncell = (counts[:, 0] > 0).sum()
+ if ncell == 0:
+ return None, None, None, None, None, None, None, 0
+ nongb_counts = np.sum(pdata[:, 1::nval], axis=0)
+ hmin = np.min(pdata[:, 2::nval], axis=0, initial=np.inf, where=counts > 0)
+ hmax = np.max(pdata[:, 3::nval], axis=0, initial=0.0, where=counts > 0)
+ hsum = np.sum(pdata[:, 4::nval], axis=0, where=counts > 0)
+ hsum2 = np.sum(pdata[:, 5::nval], axis=0, where=counts > 0)
+ counts = counts.sum(axis=0)
+ filt = counts > 0
+ hmean = hsum
+ hmean[filt] /= counts[filt]
+ hmean2 = hsum2
+ hmean2[filt] /= counts[filt]
+ hstd = np.sqrt(hmean2 - hmean ** 2)
+ bins = np.arange(hmean.shape[0])
+ return (
+ bins[filt],
+ counts[filt],
+ nongb_counts[filt],
+ hmin[filt],
+ hmax[filt],
+ hmean[filt],
+ hstd[filt],
+ ncell,
+ )
+
+
+def no_values(ax):
+ ax.text(0.5, 0.5, "No values", horizontalalignment="center", transform=ax.transAxes)
+ return
+
+
+hbins, hcounts, hno_ngb, hhmin, hhmax, hhmean, hhstd, hncell = get_total_stats(
+ hdata[:, : nval * (nbin - 1)]
+)
+htbins, htcounts, htno_ngb, hthmin, hthmax, hthmean, hthstd, htncell = get_total_stats(
+ hdata[:, nval * (nbin - 1) :]
+)
+sbins, scounts, sno_ngb, shmin, shmax, shmean, shstd, sncell = get_total_stats(
+ sdata[:, : nval * (nbin - 1)]
+)
+stbins, stcounts, stno_ngb, sthmin, sthmax, sthmean, sthstd, stncell = get_total_stats(
+ sdata[:, nval * (nbin - 1) :]
+)
+bbins, bcounts, bno_ngb, bhmin, bhmax, bhmean, bhstd, bncell = get_total_stats(
+ bdata[:, : nval * (nbin - 1)]
+)
+btbins, btcounts, btno_ngb, bthmin, bthmax, bthmean, bthstd, btncell = get_total_stats(
+ bdata[:, nval * (nbin - 1) :]
+)
+
+fig, ax = pl.subplots(3, 3, sharex=True, figsize=(10, 8))
+
+if hncell > 0:
+ ax[0][0].axhline(y=htcounts[0], linestyle="--", color="k")
+ ax[1][0].axhline(y=hthmin, linestyle=":", color="k")
+ ax[1][0].axhline(y=hthmean - hthstd, linestyle="--", color="k")
+ ax[1][0].axhline(y=hthmean, linestyle="-", color="k")
+ ax[1][0].axhline(y=hthmean + hthstd, linestyle="--", color="k")
+ ax[1][0].axhline(y=hthmax, linestyle=":", color="k")
+ ax[2][0].axhline(y=htno_ngb[0], linestyle="--", color="k")
+
+ ax[0][0].bar(hbins, hcounts, log=True)
+
+ ax[1][0].plot(hbins, hhmin, label="min")
+ ax[1][0].plot(hbins, hhmax, label="max")
+ ax[1][0].fill_between(hbins, hhmean - hhstd, hhmean + hhstd, color="C2", alpha=0.5)
+ ax[1][0].plot(hbins, hhmean, color="C2", label="mean")
+ ax[1][0].set_yscale("log")
+
+ if hno_ngb.sum() > 0:
+ ax[2][0].bar(hbins, hno_ngb, log=True)
+ else:
+ no_values(ax[2][0])
+else:
+ no_values(ax[0][0])
+ no_values(ax[1][0])
+ no_values(ax[2][0])
+
+if sncell > 0:
+ ax[0][1].axhline(y=stcounts[0], linestyle="--", color="k")
+ ax[1][1].axhline(y=sthmin, linestyle=":", color="k")
+ ax[1][1].axhline(y=sthmean - sthstd, linestyle="--", color="k")
+ ax[1][1].axhline(y=sthmean, linestyle="-", color="k")
+ ax[1][1].axhline(y=sthmean + sthstd, linestyle="--", color="k")
+ ax[1][1].axhline(y=sthmax, linestyle=":", color="k")
+ ax[2][1].axhline(y=stno_ngb[0], linestyle="--", color="k")
+
+ ax[0][1].bar(sbins, scounts, log=True)
+
+ ax[1][1].plot(sbins, shmin)
+ ax[1][1].plot(sbins, shmax)
+ ax[1][1].fill_between(sbins, shmean - shstd, shmean + shstd, color="C2", alpha=0.5)
+ ax[1][1].plot(sbins, shmean, color="C2")
+ ax[1][1].set_yscale("log")
+
+ if sno_ngb.sum() > 0:
+ ax[2][1].bar(sbins, sno_ngb, log=True)
+ else:
+ no_values(ax[2][1])
+else:
+ no_values(ax[0][1])
+ no_values(ax[1][1])
+ no_values(ax[2][1])
+
+if bncell > 0:
+ ax[0][2].axhline(y=btcounts[0], linestyle="--", color="k")
+ ax[1][2].axhline(y=bthmin, linestyle=":", color="k")
+ ax[1][2].axhline(y=bthmean - bthstd, linestyle="--", color="k")
+ ax[1][2].axhline(y=bthmean, linestyle="-", color="k")
+ ax[1][2].axhline(y=bthmean + bthstd, linestyle="--", color="k")
+ ax[1][2].axhline(y=bthmax, linestyle=":", color="k")
+ ax[2][2].axhline(y=btno_ngb[0], linestyle="--", color="k")
+
+ ax[0][2].bar(bbins, bcounts, log=True)
+
+ ax[1][2].plot(bbins, bhmin)
+ ax[1][2].plot(bbins, bhmax)
+ ax[1][2].fill_between(bbins, bhmean - bhstd, bhmean + bhstd, color="C2", alpha=0.5)
+ ax[1][2].plot(bbins, bhmean, color="C2")
+ ax[1][2].set_yscale("log")
+
+ if bno_ngb.sum() > 0:
+ ax[2][2].bar(bbins, bno_ngb, log=True)
+ else:
+ no_values(ax[2][2])
+else:
+ no_values(ax[0][2])
+ no_values(ax[1][2])
+ no_values(ax[2][2])
+
+ax[1][0].set_xlabel("iteration")
+ax[1][1].set_xlabel("iteration")
+ax[1][2].set_xlabel("iteration")
+
+ax[0][0].set_ylabel("count")
+ax[1][0].set_ylabel("h")
+ax[2][0].set_ylabel("no ngb count")
+
+ax[1][0].legend(loc="best")
+ax[1][1].plot([], [], "k-", label="converged value")
+ax[1][1].legend(loc="best")
+
+ax[0][0].set_xlim(-0.5, nbin - 1.5)
+ax[0][0].set_xticks(np.arange(nbin))
+ax[0][0].set_title("hydro")
+ax[0][1].set_title("stars")
+ax[0][2].set_title("black holes")
+ax[1][0].set_title("{0} cells".format(hncell))
+ax[1][1].set_title("{0} cells".format(sncell))
+ax[1][2].set_title("{0} cells".format(bncell))
+
+if not args.label is None:
+ pl.suptitle(args.label)
+else:
+ pl.suptitle(",".join(args.input))
+
+pl.tight_layout()
+pl.savefig(args.output, dpi=300)
--
GitLab