diff --git a/README b/README
index b0a21c610f7ab4c380e523b51c6f74b81f98a2f2..a2808dbf77e29bf304f320d927870a557597bb07 100644
--- a/README
+++ b/README
@@ -26,6 +26,8 @@ Valid options are:
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements
-g Run with an external gravitational potential
-G Run with self-gravity
+ -n {int} Execute a fixed number of time steps. Defaults to -1, which means
+ use the time_end parameter to stop.
-s Run with SPH
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-v [12] Increase the level of verbosity 1: MPI-rank 0 writes
diff --git a/configure.ac b/configure.ac
index 497107121753f212dd8b07f5a8e8eed7acdf82b5..23fbf5c6e4c4211379986a4a9c951892d044713c 100644
--- a/configure.ac
+++ b/configure.ac
@@ -72,9 +72,6 @@ if test "$enable_ipo" = "yes"; then
fi
fi
-# Add libtool support.
-LT_INIT
-
# Check for MPI. Need to do this before characterising the compiler (C99 mode),
# as this changes the compiler.
# We should consider using AX_PROG_CC_MPI to replace AC_PROG_CC when compiling
@@ -151,6 +148,9 @@ AM_CONDITIONAL([HAVEMPI],[test $enable_mpi = "yes"])
# Indicate that MPIRUN can be modified by an environement variable
AC_ARG_VAR(MPIRUN, Path to the mpirun command if non-standard)
+# Add libtool support (now that CC is defined).
+LT_INIT
+
# Need C99 and inline support.
AC_PROG_CC_C99
AC_C_INLINE
@@ -179,7 +179,7 @@ if test "$enable_opt" = "yes" ; then
ac_test_CFLAGS="yes"
CFLAGS="$old_CFLAGS $CFLAGS"
- # Check SSE & AVX support (some overlap with AX_CC_MAXOPT).
+ # Check SSE & AVX support (some overlap with AX_CC_MAXOPT).
# Don't use the SIMD_FLAGS result with Intel compilers. The -x<code>
# value from AX_CC_MAXOPT should be sufficient.
AX_EXT
@@ -287,12 +287,75 @@ AC_SUBST([METIS_LIBS])
AC_SUBST([METIS_INCS])
AM_CONDITIONAL([HAVEMETIS],[test -n "$METIS_LIBS"])
-# # Check for zlib.
-# AC_CHECK_LIB([z],[gzopen],[
-# AC_DEFINE([HAVE_LIBZ],[1],[Set to 1 if zlib is installed.])
-# LDFLAGS="$LDFLAGS -lz"
-# ],[])
+# Check for tcmalloc a fast malloc that is part of the gperftools.
+have_tcmalloc="no"
+AC_ARG_WITH([tcmalloc],
+ [AS_HELP_STRING([--with-tcmalloc],
+ [use tcmalloc library or specify the directory with lib @<:@yes/no@:>@]
+ )],
+ [with_tcmalloc="$withval"],
+ [with_tcmalloc="no"]
+)
+if test "x$with_tcmalloc" != "xno"; then
+ if test "x$with_tcmalloc" != "xyes" && test "x$with_tcmalloc" != "x"; then
+ tclibs="-L$with_tcmalloc -ltcmalloc"
+ else
+ tclibs="-ltcmalloc"
+ fi
+ AC_CHECK_LIB([tcmalloc],[tc_cfree],[have_tcmalloc="yes"],[have_tcmalloc="no"],
+ $tclibs)
+
+ # Could just have the minimal version.
+ if test "$have_tcmalloc" = "no"; then
+ if test "x$with_tcmalloc" != "xyes" && test "x$with_tcmalloc" != "x"; then
+ tclibs="-L$with_tcmalloc -ltcmalloc_minimal"
+ else
+ tclibs="-ltcmalloc_minimal"
+ fi
+ AC_CHECK_LIB([tcmalloc],[tc_cfree],[have_tcmalloc="yes"],[have_tcmalloc="no"],
+ $tclibs)
+ fi
+ if test "$have_tcmalloc" = "yes"; then
+ TCMALLOC_LIBS="$tclibs"
+
+ # These are recommended for GCC.
+ if test "$ax_cv_c_compiler_vendor" = "gnu"; then
+ CFLAGS="$CFLAGS -fno-builtin-malloc -fno-builtin-calloc -fno-builtin-realloc -fno-builtin-free"
+ fi
+ else
+ TCMALLOC_LIBS=""
+ fi
+fi
+AC_SUBST([TCMALLOC_LIBS])
+AM_CONDITIONAL([HAVETCMALLOC],[test -n "$TCMALLOC_LIBS"])
+
+# Check for -lprofiler usually part of the gpreftools along with tcmalloc.
+have_profiler="no"
+AC_ARG_WITH([profiler],
+ [AS_HELP_STRING([--with-profiler],
+ [use cpu profiler library or specify the directory with lib @<:@yes/no@:>@]
+ )],
+ [with_profiler="$withval"],
+ [with_profiler="yes"]
+)
+if test "x$with_profiler" != "xno"; then
+ if test "x$with_profiler" != "xyes" && test "x$with_profiler" != "x"; then
+ proflibs="-L$with_profiler -lprofiler"
+ else
+ proflibs="-lprofiler"
+ fi
+ AC_CHECK_LIB([profiler],[ProfilerFlush],[have_profiler="yes"],[have_profiler="no"],
+ $proflibs)
+
+ if test "$have_profiler" = "yes"; then
+ PROFILER_LIBS="$proflibs"
+ else
+ PROFILER_LIBS=""
+ fi
+fi
+AC_SUBST([PROFILER_LIBS])
+AM_CONDITIONAL([HAVEPROFILER],[test -n "$PROFILER_LIBS"])
# Check for HDF5. This is required.
AX_LIB_HDF5
@@ -410,6 +473,8 @@ AC_MSG_RESULT([
- parallel : $have_parallel_hdf5
Metis enabled : $have_metis
libNUMA enabled : $have_numa
+ Using tcmalloc : $have_tcmalloc
+ CPU profiler : $have_profiler
])
# Generate output.
diff --git a/examples/Makefile.am b/examples/Makefile.am
index 735817c24cc52786e4c562e46e3619fb4a9a2e34..22911e903e8bbf8265d712185f41de2bf4f53553 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -21,13 +21,17 @@ MYFLAGS = -DTIMER
# Add the source directory and debug to CFLAGS
AM_CFLAGS = -I../src $(HDF5_CPPFLAGS)
-AM_LDFLAGS =
+AM_LDFLAGS = $(HDF5_LDFLAGS)
-MPI_THREAD_LIBS = @MPI_THREAD_LIBS@
+# Extra libraries.
+EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS)
+
+# MPI libraries.
MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
-# Set-up the library
+
+# Programs.
bin_PROGRAMS = swift swift_fixdt
# Build MPI versions as well?
@@ -45,20 +49,20 @@ endif
# Sources for swift
swift_SOURCES = main.c
swift_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_keep $(ENGINE_POLICY_SETAFFINITY)"
-swift_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+swift_LDADD = ../src/.libs/libswiftsim.a $(EXTRA_LIBS)
swift_fixdt_SOURCES = main.c
swift_fixdt_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep $(ENGINE_POLICY_SETAFFINITY)"
-swift_fixdt_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
+swift_fixdt_LDADD = ../src/.libs/libswiftsim.a $(EXTRA_LIBS)
# Sources for swift_mpi, do we need an affinity policy for MPI?
swift_mpi_SOURCES = main.c
swift_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_keep $(ENGINE_POLICY_SETAFFINITY)"
-swift_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
+swift_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(MPI_LIBS) $(EXTRA_LIBS)
swift_fixdt_mpi_SOURCES = main.c
swift_fixdt_mpi_CFLAGS = $(MYFLAGS) $(AM_CFLAGS) $(MPI_FLAGS) -DENGINE_POLICY="engine_policy_fixdt | engine_policy_keep $(ENGINE_POLICY_SETAFFINITY)"
-swift_fixdt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(MPI_LIBS)
+swift_fixdt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(MPI_LIBS) $(EXTRA_LIBS)
# Scripts to generate ICs
EXTRA_DIST = UniformBox/makeIC.py UniformBox/run.sh UniformBox/uniformBox.yml \
diff --git a/examples/main.c b/examples/main.c
index fa0ac2567ac253568903ce8c47b2c9c38fee3704..9907504e5032e2ba72ef3249dc66149b2a5eeefa 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -74,6 +74,7 @@ void print_help_message() {
printf(" %2s %8s %s\n", "-g", "",
"Run with an external gravitational potential");
printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity");
+ printf(" %2s %8s %s\n", "-n", "{int}", "Execute a fixed number of time steps");
printf(" %2s %8s %s\n", "-s", "", "Run with SPH");
printf(" %2s %8s %s\n", "-t", "{int}",
"The number of threads to use on each MPI rank. Defaults to 1 if not "
@@ -138,6 +139,7 @@ int main(int argc, char *argv[]) {
int with_aff = 0;
int dry_run = 0;
int dump_tasks = 0;
+ int nsteps = -1;
int with_cosmology = 0;
int with_external_gravity = 0;
int with_self_gravity = 0;
@@ -150,7 +152,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
- while ((c = getopt(argc, argv, "acdef:gGhst:v:y:")) != -1) switch (c) {
+ while ((c = getopt(argc, argv, "acdef:gGhn:st:v:y:")) != -1) switch (c) {
case 'a':
with_aff = 1;
break;
@@ -179,6 +181,13 @@ int main(int argc, char *argv[]) {
case 'h':
if (myrank == 0) print_help_message();
return 0;
+ case 'n':
+ if (sscanf(optarg, "%d", &nsteps) != 1) {
+ if (myrank == 0) printf("Error parsing fixed number of steps.\n");
+ if (myrank == 0) print_help_message();
+ return 1;
+ }
+ break;
case 's':
with_hydro = 1;
break;
@@ -322,18 +331,21 @@ int main(int argc, char *argv[]) {
size_t Ngas = 0, Ngpart = 0;
double dim[3] = {0., 0., 0.};
int periodic = 0;
+ int flag_entropy_ICs = 0;
if (myrank == 0) clocks_gettime(&tic);
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
- myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+ &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL, dry_run);
#else
read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
- myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+ &flag_entropy_ICs, myrank, nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL, dry_run);
#endif
#else
read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
- dry_run);
+ &flag_entropy_ICs, dry_run);
#endif
if (myrank == 0) {
clocks_gettime(&toc);
@@ -354,7 +366,7 @@ int main(int argc, char *argv[]) {
free(parts);
parts = NULL;
for (size_t k = 0; k < Ngpart; ++k)
- if (gparts[k].id > 0) error("Linking problem");
+ if (gparts[k].id_or_neg_offset < 0) error("Linking problem");
Ngas = 0;
}
@@ -466,7 +478,7 @@ int main(int argc, char *argv[]) {
#endif
/* Initialise the particles */
- engine_init_particles(&e);
+ engine_init_particles(&e, flag_entropy_ICs);
/* Legend */
if (myrank == 0)
@@ -474,7 +486,7 @@ int main(int argc, char *argv[]) {
"Updates", "g-Updates", "Wall-clock time", clocks_getunit());
/* Main simulation loop */
- for (int j = 0; !engine_is_done(&e); j++) {
+ for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
/* Repartition the space amongst the nodes? */
#ifdef WITH_MPI
diff --git a/src/Makefile.am b/src/Makefile.am
index 21b5ef25ca21202caaa9107bfbf09e62aa66011d..bc263d6fb0b2a8378e26a3627cf26204c7b09bb4 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -24,6 +24,13 @@ AM_LDFLAGS = $(HDF5_LDFLAGS) -version-info 0:0:0
# The git command, if available.
GIT_CMD = @GIT_CMD@
+# Additional dependencies for shared libraries.
+EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS)
+
+# MPI libraries.
+MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
+MPI_FLAGS = -DWITH_MPI $(METIS_INCS)
+
# Build the libswiftsim library
lib_LTLIBRARIES = libswiftsim.la
# Build a MPI-enabled version too?
@@ -65,11 +72,13 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
# Sources and flags for regular library
libswiftsim_la_SOURCES = $(AM_SOURCES)
+libswiftsim_la_CFLAGS = $(AM_CFLAGS)
+libswiftsim_la_LDFLAGS = $(AM_LDFLAGS) $(EXTRA_LIBS)
# Sources and flags for MPI library
libswiftsim_mpi_la_SOURCES = $(AM_SOURCES)
-libswiftsim_mpi_la_CFLAGS = $(AM_CFLAGS) -DWITH_MPI $(METIS_INCS)
-libswiftsim_mpi_la_LDFLAGS = $(AM_LDFLAGS) -DWITH_MPI $(METIS_LIBS)
+libswiftsim_mpi_la_CFLAGS = $(AM_CFLAGS) $(MPI_FLAGS)
+libswiftsim_mpi_la_LDFLAGS = $(AM_LDFLAGS) $(MPI_LIBS) $(EXTRA_LIBS)
libswiftsim_mpi_la_SHORTNAME = mpi
diff --git a/src/cell.c b/src/cell.c
index c9e35c1fb9beb8d011d67d0f7107f64b0143f44f..cd83ddab570a3c848e0502629534743e8f383a4f 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -222,7 +222,7 @@ int cell_pack_ti_ends(struct cell *c, int *ti_ends) {
/* Pack this cell's data. */
ti_ends[0] = c->ti_end_min;
-
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
@@ -238,7 +238,7 @@ int cell_unpack_ti_ends(struct cell *c, int *ti_ends) {
/* Unpack this cell's data. */
c->ti_end_min = ti_ends[0];
-
+
/* Fill in the progeny, depth-first recursion. */
int count = 1;
for (int k = 0; k < 8; k++)
@@ -410,9 +410,11 @@ void cell_gunlocktree(struct cell *c) {
* @brief Sort the parts into eight bins along the given pivots.
*
* @param c The #cell array to be sorted.
+ * @param parts_offset Offset of the cell parts array relative to the
+ * space's parts array, i.e. c->parts - s->parts.
*/
-void cell_split(struct cell *c) {
+void cell_split(struct cell *c, ptrdiff_t parts_offset) {
int i, j;
const int count = c->count, gcount = c->gcount;
@@ -530,8 +532,7 @@ void cell_split(struct cell *c) {
}
/* Re-link the gparts. */
- for (int k = 0; k < count; k++)
- if (parts[k].gpart != NULL) parts[k].gpart->part = &parts[k];
+ part_relink_gparts(parts, count, parts_offset);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that _all_ the parts have been assigned to a cell. */
@@ -626,8 +627,7 @@ void cell_split(struct cell *c) {
}
/* Re-link the parts. */
- for (int k = 0; k < gcount; k++)
- if (gparts[k].id > 0) gparts[k].part->gpart = &gparts[k];
+ part_relink_parts(gparts, gcount, parts - parts_offset);
}
/**
diff --git a/src/cell.h b/src/cell.h
index 24b679caf480c741de42001ff9abf7dcfc172de4..b5cc49d05c1cee6e0a884cdb708abbaff4281829 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -24,6 +24,9 @@
#define SWIFT_CELL_H
/* Includes. */
+#include <stddef.h>
+
+/* Local includes. */
#include "lock.h"
#include "multipole.h"
#include "part.h"
@@ -178,7 +181,7 @@ struct cell {
((int)(k) + (cdim)[2] * ((int)(j) + (cdim)[1] * (int)(i)))
/* Function prototypes. */
-void cell_split(struct cell *c);
+void cell_split(struct cell *c, ptrdiff_t parts_offset);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
diff --git a/src/common_io.c b/src/common_io.c
index 03d5618cd0cafe5ff6067cf72107baacffb5d85a..971fe6b01c682c489d3444ec1d47d6a902250bb8 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -516,12 +516,10 @@ void prepare_dm_gparts(struct gpart* const gparts, size_t Ndm) {
/* Let's give all these gparts a negative id */
for (size_t i = 0; i < Ndm; ++i) {
-
/* 0 or negative ids are not allowed */
- if (gparts[i].id <= 0)
- error("0 or negative ID for DM particle %zd: ID=%lld", i, gparts[i].id);
-
- gparts[i].id = -gparts[i].id;
+ if (gparts[i].id_or_neg_offset <= 0)
+ error("0 or negative ID for DM particle %zd: ID=%lld", i,
+ gparts[i].id_or_neg_offset);
}
}
@@ -555,7 +553,7 @@ void duplicate_hydro_gparts(struct part* const parts,
gparts[i + Ndm].mass = parts[i].mass;
/* Link the particles */
- gparts[i + Ndm].part = &parts[i];
+ gparts[i + Ndm].id_or_neg_offset = -i;
parts[i].gpart = &gparts[i + Ndm];
}
}
@@ -580,9 +578,8 @@ void collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
* gparts[i].part); */
/* And collect the DM ones */
- if (gparts[i].id < 0LL) {
- memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
- dmparts[count].id = -dmparts[count].id;
+ if (gparts[i].id_or_neg_offset > 0) {
+ dmparts[count] = gparts[i];
count++;
}
}
diff --git a/src/debug.c b/src/debug.c
index 77238ab135e3f27b8c2c43b0ec18e7745493b138..e7b60c52f92fdb64fdefc8ac5c5e353f1ab57e6e 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -59,8 +59,8 @@
*
* (Should be used for debugging only as it runs in O(N).)
*/
-void printParticle(struct part *parts, struct xpart *xparts, long long int id,
- size_t N) {
+void printParticle(const struct part *parts, struct xpart *xparts,
+ long long int id, size_t N) {
int found = 0;
@@ -82,24 +82,27 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
* the standard output.
*
* @param gparts The array of g-particles.
+ * @param parts The array of particles.
* @param id The id too look for.
* @param N The size of the array of g-particles.
*
* (Should be used for debugging only as it runs in O(N).)
*/
-void printgParticle(struct gpart *gparts, long long int id, size_t N) {
+void printgParticle(const struct gpart *gparts, const struct part *parts,
+ long long int id, size_t N) {
int found = 0;
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
- if (gparts[i].id == -id) {
- printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
+ if (gparts[i].id_or_neg_offset == id) {
+ printf("## gParticle[%zd] (DM) :\n id=%lld", i, id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
- } else if (gparts[i].id > 0 && gparts[i].part->id == id) {
- printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
+ } else if (gparts[i].id_or_neg_offset < 0 &&
+ parts[-gparts[i].id_or_neg_offset].id == id) {
+ printf("## gParticle[%zd] (hydro) :\n id=%lld", i, id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
@@ -114,9 +117,9 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
* @param p The particle to print
* @param xp The extended data ot the particle to print
*/
-void printParticle_single(struct part *p, struct xpart *xp) {
+void printParticle_single(const struct part *p, const struct xpart *xp) {
- printf("## Particle: id=%lld ", p->id);
+ printf("## Particle: id=%lld", p->id);
hydro_debug_particle(p, xp);
printf("\n");
}
@@ -128,7 +131,7 @@ void printParticle_single(struct part *p, struct xpart *xp) {
*/
void printgParticle_single(struct gpart *gp) {
- printf("## g-Particle: id=%lld ", gp->id);
+ printf("## g-Particle: id=%lld ", gp->id_or_neg_offset);
gravity_debug_particle(gp);
printf("\n");
}
diff --git a/src/debug.h b/src/debug.h
index 13be15adb867e8bafe95e2900f68caaa36192510..367241201977d9b79a8c2913dbae5d08f1148529 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -19,14 +19,15 @@
#ifndef SWIFT_DEBUG_H
#define SWIFT_DEBUG_H
-struct part;
-struct gpart;
-struct xpart;
+/* Includes. */
+#include "cell.h"
+#include "part.h"
-void printParticle(struct part *parts, struct xpart *xparts, long long int id,
- size_t N);
-void printgParticle(struct gpart *parts, long long int id, size_t N);
-void printParticle_single(struct part *p, struct xpart *xp);
+void printParticle(const struct part *parts, struct xpart *xparts,
+ long long int id, size_t N);
+void printgParticle(const struct gpart *gparts, const struct part *parts,
+ long long int id, size_t N);
+void printParticle_single(const struct part *p, const struct xpart *xp);
void printgParticle_single(struct gpart *gp);
#ifdef HAVE_METIS
diff --git a/src/engine.c b/src/engine.c
index cd8eb6d3c8f24d75f43635ed590480b7686c2434..aa78112329aca492eb77f256ec0ccd939bcbda32 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -102,8 +102,8 @@ struct link *engine_addlink(struct engine *e, struct link *l, struct task *t) {
}
/**
- * @brief Generate the hierarchical tasks for a hierarchy of cells - i.e. all
- * the O(Npart) tasks.
+ * @brief Generate the gravity hierarchical tasks for a hierarchy of cells -
+ * i.e. all the O(Npart) tasks.
*
* Tasks are only created here. The dependencies will be added later on.
*
@@ -111,8 +111,8 @@ struct link *engine_addlink(struct engine *e, struct link *l, struct task *t) {
* @param c The #cell.
* @param super The super #cell.
*/
-void engine_make_hierarchical_tasks(struct engine *e, struct cell *c,
- struct cell *super) {
+void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c,
+ struct cell *super) {
struct scheduler *s = &e->sched;
const int is_with_external_gravity =
@@ -132,36 +132,91 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c,
if (c->nodeID == e->nodeID) {
/* Add the init task. */
- c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0, c,
- NULL, 0);
+ if (c->init == NULL)
+ c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0,
+ c, NULL, 0);
/* Add the drift task. */
- c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
- c, NULL, 0);
+ if (c->drift == NULL)
+ c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0,
+ 0, c, NULL, 0);
/* Add the kick task that matches the policy. */
if (is_fixdt) {
- c->kick = scheduler_addtask(s, task_type_kick_fixdt, task_subtype_none,
- 0, 0, c, NULL, 0);
+ if (c->kick == NULL)
+ c->kick = scheduler_addtask(s, task_type_kick_fixdt,
+ task_subtype_none, 0, 0, c, NULL, 0);
} else {
- c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0,
- c, NULL, 0);
+ if (c->kick == NULL)
+ c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0,
+ 0, c, NULL, 0);
}
- if (c->count > 0) {
+ if (is_with_external_gravity)
+ c->grav_external = scheduler_addtask(
+ s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
+ }
+ }
- /* Generate the ghost task. */
- c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
- 0, c, NULL, 0);
- }
+ /* Set the super-cell. */
+ c->super = super;
+
+ /* Recurse. */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL)
+ engine_make_gravity_hierarchical_tasks(e, c->progeny[k], super);
+}
- if (c->gcount > 0) {
+/**
+ * @brief Generate the hydro hierarchical tasks for a hierarchy of cells -
+ * i.e. all the O(Npart) tasks.
+ *
+ * Tasks are only created here. The dependencies will be added later on.
+ *
+ * @param e The #engine.
+ * @param c The #cell.
+ * @param super The super #cell.
+ */
+void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
+ struct cell *super) {
+
+ struct scheduler *s = &e->sched;
+ const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
- /* Add the external gravity tasks */
- if (is_with_external_gravity)
- c->grav_external = scheduler_addtask(
- s, task_type_grav_external, task_subtype_none, 0, 0, c, NULL, 0);
+ /* Is this the super-cell? */
+ if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
+
+ /* This is the super cell, i.e. the first with density tasks attached. */
+ super = c;
+
+ /* Local tasks only... */
+ if (c->nodeID == e->nodeID) {
+
+ /* Add the init task. */
+ if (c->init == NULL)
+ c->init = scheduler_addtask(s, task_type_init, task_subtype_none, 0, 0,
+ c, NULL, 0);
+
+ /* Add the drift task. */
+ if (c->drift == NULL)
+ c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0,
+ 0, c, NULL, 0);
+
+ /* Add the kick task that matches the policy. */
+ if (is_fixdt) {
+ if (c->kick == NULL)
+ c->kick = scheduler_addtask(s, task_type_kick_fixdt,
+ task_subtype_none, 0, 0, c, NULL, 0);
+ } else {
+ if (c->kick == NULL)
+ c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0,
+ 0, c, NULL, 0);
}
+
+ /* Generate the ghost task. */
+ c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
+ c, NULL, 0);
}
}
@@ -172,7 +227,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c,
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
- engine_make_hierarchical_tasks(e, c->progeny[k], super);
+ engine_make_hydro_hierarchical_tasks(e, c->progeny[k], super);
}
/**
@@ -270,11 +325,11 @@ void engine_redistribute(struct engine *e) {
}
#ifdef SWIFT_DEBUG_CHECKS
- if (s->parts[k].gpart->id < 0)
+ if (s->parts[k].gpart->id_or_neg_offset >= 0)
error("Trying to link a partnerless gpart !");
#endif
- s->parts[k].gpart->id = count_this_dest;
+ s->parts[k].gpart->id_or_neg_offset = -count_this_dest;
count_this_dest++;
}
}
@@ -458,13 +513,14 @@ void engine_redistribute(struct engine *e) {
for (size_t k = offset_gparts; k < offset_gparts + count_gparts; ++k) {
/* Does this gpart have a partner ? */
- if (gparts_new[k].id >= 0) {
+ if (gparts_new[k].id_or_neg_offset <= 0) {
- const size_t partner_index = offset_parts + gparts_new[k].id;
+ const ptrdiff_t partner_index =
+ offset_parts - gparts_new[k].id_or_neg_offset;
/* Re-link */
- gparts_new[k].part = &parts_new[partner_index];
- gparts_new[k].part->gpart = &gparts_new[k];
+ gparts_new[k].id_or_neg_offset = -partner_index;
+ parts_new[partner_index].gpart = &gparts_new[k];
}
}
@@ -485,22 +541,22 @@ void engine_redistribute(struct engine *e) {
/* Verify that the links are correct */
for (size_t k = 0; k < nr_gparts; ++k) {
- if (gparts_new[k].id > 0) {
+ if (gparts_new[k].id_or_neg_offset <= 0) {
- if (gparts_new[k].part->gpart != &gparts_new[k])
- error("Linking problem !");
+ struct part *part = &parts_new[-gparts_new[k].id_or_neg_offset];
- if (gparts_new[k].x[0] != gparts_new[k].part->x[0] ||
- gparts_new[k].x[1] != gparts_new[k].part->x[1] ||
- gparts_new[k].x[2] != gparts_new[k].part->x[2])
+ if (part->gpart != &gparts_new[k]) error("Linking problem !");
+
+ if (gparts_new[k].x[0] != part->x[0] ||
+ gparts_new[k].x[1] != part->x[1] || gparts_new[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < nr_parts; ++k) {
- if (parts_new[k].gpart != NULL) {
-
- if (parts_new[k].gpart->part != &parts_new[k]) error("Linking problem !");
+ if (parts_new[k].gpart != NULL &&
+ parts_new[k].gpart->id_or_neg_offset != -k) {
+ error("Linking problem !");
}
}
#endif
@@ -682,7 +738,6 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj,
* @param t_xv The recv_xv #task, if it has already been created.
* @param t_rho The recv_rho #task, if it has already been created.
*/
-
void engine_addtasks_recv(struct engine *e, struct cell *c, struct task *t_xv,
struct task *t_rho, struct task *t_ti) {
@@ -748,7 +803,8 @@ void engine_exchange_cells(struct engine *e) {
MPI_Status status;
ticks tic = getticks();
- /* Run through the cells and get the size of the ones that will be sent */
+ /* Run through the cells and get the size of the ones that will be sent off.
+ */
int count_out = 0;
for (int k = 0; k < nr_cells; k++) {
offset[k] = count_out;
@@ -913,7 +969,8 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the associated gpart with the buffer offset of the part. */
if (s->parts[offset_parts + k].gpart != NULL) {
- s->parts[offset_parts + k].gpart->id = e->proxies[pid].nr_parts_out;
+ s->parts[offset_parts + k].gpart->id_or_neg_offset =
+ -e->proxies[pid].nr_parts_out;
}
/* Load the part and xpart into the proxy. */
@@ -929,7 +986,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
error(
"Do not have a proxy for the requested nodeID %i for part with "
"id=%lli, x=[%e,%e,%e].",
- node_id, s->gparts[offset_parts + k].id,
+ node_id, s->gparts[offset_parts + k].id_or_neg_offset,
s->gparts[offset_gparts + k].x[0], s->gparts[offset_parts + k].x[1],
s->gparts[offset_gparts + k].x[2]);
proxy_gparts_load(&e->proxies[pid], &s->gparts[offset_gparts + k], 1);
@@ -989,7 +1046,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
s->xparts = xparts_new;
for (size_t k = 0; k < offset_parts; k++) {
if (s->parts[k].gpart != NULL) {
- s->parts[k].gpart->part = &s->parts[k];
+ s->parts[k].gpart->id_or_neg_offset = -k;
}
}
}
@@ -1004,8 +1061,8 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
free(s->gparts);
s->gparts = gparts_new;
for (size_t k = 0; k < offset_gparts; k++) {
- if (s->gparts[k].id > 0) {
- s->gparts[k].part->gpart = &s->gparts[k];
+ if (s->gparts[k].id_or_neg_offset < 0) {
+ s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
}
}
@@ -1078,9 +1135,10 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
/* Re-link the gparts. */
for (int k = 0; k < p->nr_gparts_in; k++) {
struct gpart *gp = &s->gparts[offset_gparts + count_gparts + k];
- if (gp->id >= 0) {
- struct part *p = &s->parts[offset_parts + count_parts + gp->id];
- gp->part = p;
+ if (gp->id_or_neg_offset <= 0) {
+ struct part *p =
+ &s->parts[offset_gparts + count_parts - gp->id_or_neg_offset];
+ gp->id_or_neg_offset = s->parts - p;
p->gpart = gp;
}
}
@@ -1504,8 +1562,14 @@ void engine_maketasks(struct engine *e) {
engine_count_and_link_tasks(e);
/* Append hierarchical tasks to each cells */
- for (int k = 0; k < nr_cells; k++)
- engine_make_hierarchical_tasks(e, &cells[k], NULL);
+ if (e->policy & engine_policy_hydro)
+ for (int k = 0; k < nr_cells; k++)
+ engine_make_hydro_hierarchical_tasks(e, &cells[k], NULL);
+
+ if ((e->policy & engine_policy_self_gravity) ||
+ (e->policy & engine_policy_external_gravity))
+ for (int k = 0; k < nr_cells; k++)
+ engine_make_gravity_hierarchical_tasks(e, &cells[k], NULL);
/* Run through the tasks and make force tasks for each density task.
Each force task depends on the cell ghosts and unlocks the kick task
@@ -1805,7 +1869,6 @@ void engine_print_task_counts(struct engine *e) {
*
* @param e The #engine.
*/
-
void engine_rebuild(struct engine *e) {
const ticks tic = getticks();
@@ -2183,8 +2246,10 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
*forward in time.
*
* @param e The #engine
+ * @param flag_entropy_ICs Did the 'Internal Energy' of the particles actually
+ *contain entropy ?
*/
-void engine_init_particles(struct engine *e) {
+void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
struct space *s = e->s;
@@ -2252,7 +2317,7 @@ void engine_init_particles(struct engine *e) {
TIMER_TOC(timer_runners);
/* Apply some conversions (e.g. internal energy -> entropy) */
- space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
+ if (!flag_entropy_ICs) space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
clocks_gettime(&time2);
@@ -2537,8 +2602,7 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
s->xparts = xparts_new;
/* Re-link the gparts. */
- for (size_t k = 0; k < s->nr_parts; k++)
- if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
+ part_relink_gparts(s->parts, s->nr_parts, 0);
/* Re-allocate the local gparts. */
if (e->verbose)
@@ -2554,31 +2618,28 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
s->gparts = gparts_new;
/* Re-link the parts. */
- for (size_t k = 0; k < s->nr_gparts; k++)
- if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
+ part_relink_parts(s->gparts, s->nr_gparts, s->parts);
#ifdef SWIFT_DEBUG_CHECKS
-
/* Verify that the links are correct */
for (size_t k = 0; k < s->nr_gparts; ++k) {
- if (s->gparts[k].id > 0) {
+ if (s->gparts[k].id_or_neg_offset <= 0) {
- if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+ struct part *part = &s->parts[-s->gparts[k].id_or_neg_offset];
- if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
- s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
- s->gparts[k].x[2] != s->gparts[k].part->x[2])
+ if (part->gpart != &s->gparts[k]) error("Linking problem !");
+
+ if (s->gparts[k].x[0] != part->x[0] || s->gparts[k].x[1] != part->x[1] ||
+ s->gparts[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < s->nr_parts; ++k) {
- if (s->parts[k].gpart != NULL) {
- if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
- }
+ if (s->parts[k].gpart != NULL && s->parts[k].gpart->id_or_neg_offset != -k)
+ error("Linking problem !");
}
-
#endif
#else
diff --git a/src/engine.h b/src/engine.h
index d1bfecc568355a9cf5aa57590ccba1df81b05b8a..3dfe22845a66ba02fcc270e13647efd4390c7256 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -214,7 +214,7 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e);
void engine_print(struct engine *e);
-void engine_init_particles(struct engine *e);
+void engine_init_particles(struct engine *e, int flag_entropy_ICs);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, struct partition *initial_partition);
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 21d775d703c8bd7000ad03386fe818124642a8f1..7cf375a1fdf7bccc4131dc415ab2d4acbbf2d3bc 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void gravity_debug_particle(
- struct gpart* p) {
+ const struct gpart* p) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 129c4b39828ca73d2d80d79edbdaa8ec4d5a9e01..74f364dd97361f0513755bedec83fe7cb277c36b 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -39,8 +39,8 @@ __attribute__((always_inline)) INLINE static void darkmatter_read_particles(
COMPULSORY);
readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v_full,
COMPULSORY);
- readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset, id,
- COMPULSORY);
+ readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset,
+ id_or_neg_offset, COMPULSORY);
}
/**
@@ -75,6 +75,6 @@ __attribute__((always_inline)) INLINE static void darkmatter_write_particles(
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, v_full, us,
UNIT_CONV_SPEED);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
- ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset, id, us,
- UNIT_CONV_NO_UNITS);
+ ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset,
+ id_or_neg_offset, us, UNIT_CONV_NO_UNITS);
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index a54fc1d18e7e84a1ca7a7e7913117ed207549271..b6c9e62559207cb25323c8a108e4bffc87ea0fcf 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -47,14 +47,8 @@ struct gpart {
/* float tx; */
/* float tv; */
- /* Anonymous union for id/part. */
- union {
-
- /* Particle ID. */
- long long id;
-
- /* Pointer to corresponding SPH part. */
- struct part* part;
- };
+ /* Particle ID. If negative, it is the negative offset of the #part with
+ which this gpart is linked. */
+ long long id_or_neg_offset;
} __attribute__((aligned(gpart_align)));
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 79ee392d46ca75a2c097bf045b2d82c9f3dc96c0..9b8136e8349398e4e9e804459a5de23d21043564 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
- struct part* p, struct xpart* xp) {
+ const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index de71963b05149fa0b2df6222424478a6ad9b1f44..fb1f7af922b2272779c7251b15e304880a2af520 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -122,3 +122,10 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
const_max_u_change);
}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 31e89d438bc96c0f0f2ba56249664d28036cb607..b67e79182ccaaee7c0421c57a91ec9fa2adae65c 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
- struct part* p, struct xpart* xp) {
+ const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index b977f25386fab0787c925635297a48fa85a8df24..f7ee04cbbad07ebb06d00ee39a614b9fed5c4dfd 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -83,8 +83,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "SmoothingLength",
FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, us,
UNIT_CONV_LENGTH);
- writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy", FLOAT, N,
- 1, parts, N_total, mpi_rank, offset, entropy, us,
+ writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "InternalEnergy",
+ FLOAT, N, 1, parts, N_total, mpi_rank, offset, entropy, us,
UNIT_CONV_ENTROPY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id, us,
@@ -113,3 +113,10 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 1; }
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 127ba75e99418b6a5dc197a44ccdc77de3cdef15..1e0e3f1f111149a3babeef431893dff2cbe6bb3c 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -18,7 +18,7 @@
******************************************************************************/
__attribute__((always_inline)) INLINE static void hydro_debug_particle(
- struct part* p, struct xpart* xp) {
+ const struct part* p, const struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 1bbfe1065358017e0b27e2131ff717848221fe9c..bed6b8f8221dc3a6e94c65d101c68aa9b3bdea91 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -115,3 +115,10 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
const_max_u_change);
}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
diff --git a/src/parallel_io.c b/src/parallel_io.c
index c5cac1cb5efc6e533e599867e39cdd7c7b2c87fa..66b02d2bff2bc0d6f2dd69d994a93576e8b96583 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -367,8 +367,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
*/
void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
- int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info, int dry_run) {
+ int* periodic, int* flag_entropy, int mpi_rank,
+ int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
@@ -404,6 +404,7 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
if (h_grp < 0) error("Error while opening file header\n");
/* Read the relevant information and print status */
+ readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -626,6 +627,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
double MassTable[6] = {0., 0., 0., 0., 0., 0.};
writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+ flagEntropy[0] = writeEntropyFlag();
writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
NUM_PARTICLE_TYPES);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 26757cb679e475d6acd2ce3c408135dfe5e49081..709479fab029e862caddb36c4216b85077e96cec 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -36,8 +36,8 @@
void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
- int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info, int dry_run);
+ int* periodic, int* flag_entropy, int mpi_rank,
+ int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run);
void write_output_parallel(struct engine* e, const char* baseName,
struct UnitSystem* us, int mpi_rank, int mpi_size,
diff --git a/src/part.c b/src/part.c
index b89abdde40fe8c7a57d1e9ac9e18fece83ba1f21..b00eaccaae0e86f7c4e8019a307f0bf455687b7c 100644
--- a/src/part.c
+++ b/src/part.c
@@ -29,6 +29,36 @@
#include "error.h"
#include "part.h"
+/**
+ * @brief Re-link the #gpart%s associated with the list of #part%s.
+ *
+ * @param parts The list of #part.
+ * @param N The number of particles to re-link;
+ * @param offset The offset of #part%s relative to the global parts list.
+ */
+void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset) {
+ for (size_t k = 0; k < N; k++) {
+ if (parts[k].gpart) {
+ parts[k].gpart->id_or_neg_offset = -(k + offset);
+ }
+ }
+}
+
+/**
+ * @brief Re-link the #gpart%s associated with the list of #part%s.
+ *
+ * @param gparts The list of #gpart.
+ * @param N The number of particles to re-link;
+ * @param parts The global part array in which to find the #gpart offsets.
+ */
+void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts) {
+ for (size_t k = 0; k < N; k++) {
+ if (gparts[k].id_or_neg_offset <= 0) {
+ parts[-gparts[k].id_or_neg_offset].gpart = &gparts[k];
+ }
+ }
+}
+
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
MPI_Datatype part_mpi_type;
diff --git a/src/part.h b/src/part.h
index e99be6e51a9bd74dd9eec8f590e80989e83ec2e1..efca7b6b5bef49f20df1e2c45b30f65ecbbf4960 100644
--- a/src/part.h
+++ b/src/part.h
@@ -22,6 +22,9 @@
/* Config parameters. */
#include "../config.h"
+/* Standard headers. */
+#include <stddef.h>
+
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
@@ -49,6 +52,8 @@
/* Import the right gravity particle definition */
#include "./gravity/Default/gravity_part.h"
+void part_relink_gparts(struct part *parts, size_t N, ptrdiff_t offset);
+void part_relink_parts(struct gpart *gparts, size_t N, struct part *parts);
#ifdef WITH_MPI
/* MPI data type for the particle transfers */
extern MPI_Datatype part_mpi_type;
diff --git a/src/runner.c b/src/runner.c
index 249479546bf98ccd30d2ac5ee484c727424fe153..349ae7cc1eb0a955127c0d85b101b085bc6264f0 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -759,7 +759,7 @@ void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart */
- if (gp->id < 0) {
+ if (gp->id_or_neg_offset > 0) {
/* First, finish the force calculation */
gravity_end_force(gp);
@@ -871,7 +871,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
- if (gp->id < 0) {
+ if (gp->id_or_neg_offset > 0) {
if (gp->ti_end <= ti_current) {
diff --git a/src/serial_io.c b/src/serial_io.c
index 7e78276dc83430655b4ea4de2fb7425e71e07966..03d3c38850e0a7800c83d7dd0942025632bf1d0b 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -424,8 +424,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
*/
void read_ic_serial(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
- int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info, int dry_run) {
+ int* periodic, int* flag_entropy, int mpi_rank,
+ int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
@@ -462,6 +462,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
if (h_grp < 0) error("Error while opening file header\n");
/* Read the relevant information and print status */
+ readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -488,6 +489,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
}
/* Now need to broadcast that information to all ranks. */
+ MPI_Bcast(flag_entropy, 1, MPI_INT, 0, comm);
MPI_Bcast(periodic, 1, MPI_INT, 0, comm);
MPI_Bcast(&N_total, NUM_PARTICLE_TYPES, MPI_LONG_LONG, 0, comm);
MPI_Bcast(dim, 3, MPI_DOUBLE, 0, comm);
@@ -699,6 +701,7 @@ void write_output_serial(struct engine* e, const char* baseName,
double MassTable[6] = {0., 0., 0., 0., 0., 0.};
writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+ flagEntropy[0] = writeEntropyFlag();
writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
NUM_PARTICLE_TYPES);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
diff --git a/src/serial_io.h b/src/serial_io.h
index 6b64624772214494a43316d3a8aa3910c7238dc8..fa88d5cb87edfba2c101b6096b8a57b04b7178cd 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -36,8 +36,8 @@
void read_ic_serial(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
- int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
- MPI_Info info, int dry_run);
+ int* periodic, int* flag_entropy, int mpi_rank,
+ int mpi_size, MPI_Comm comm, MPI_Info info, int dry_run);
void write_output_serial(struct engine* e, const char* baseName,
struct UnitSystem* us, int mpi_rank, int mpi_size,
diff --git a/src/single_io.c b/src/single_io.c
index 3f65aae0b5d495670f2b4862e466ec849f997d63..3893f1c86784e07e5cd78afc32a49339e60840ba 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -322,6 +322,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
* @param Ngas (output) number of Gas particles read.
* @param Ngparts (output) The number of #gpart read.
* @param periodic (output) 1 if the volume is periodic, 0 if not.
+ * @param flag_entropy 1 if the ICs contained Entropy in the InternalEnergy
+ * field
* @param dry_run If 1, don't read the particle. Only allocates the arrays.
*
* Opens the HDF5 file fileName and reads the particles contained
@@ -334,7 +336,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
*/
void read_ic_single(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
- int* periodic, int dry_run) {
+ int* periodic, int* flag_entropy, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
@@ -368,6 +370,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
if (h_grp < 0) error("Error while opening file header\n");
/* Read the relevant information and print status */
+ readAttribute(h_grp, "Flag_Entropy_ICs", INT, flag_entropy);
readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord);
@@ -549,6 +552,7 @@ void write_output_single(struct engine* e, const char* baseName,
double MassTable[NUM_PARTICLE_TYPES] = {0};
writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0};
+ flagEntropy[0] = writeEntropyFlag();
writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
NUM_PARTICLE_TYPES);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
diff --git a/src/single_io.h b/src/single_io.h
index d2c87655e1c91b92e8ccf2aa50d2e0bf98f13482..4936aa3eb35a5c76a5a8f242c9cc79bbc86161a0 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -31,7 +31,7 @@
void read_ic_single(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ndm,
- int* periodic, int dry_run);
+ int* periodic, int* flag_entropy, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
struct UnitSystem* us);
diff --git a/src/space.c b/src/space.c
index c4173a58895bb36ca7f81a69128574abd9700f2f..0b42a3643f5740881da85e0533c02f7d8a0ee958 100644
--- a/src/space.c
+++ b/src/space.c
@@ -442,10 +442,10 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
s->parts[k] = s->parts[nr_parts];
s->parts[nr_parts] = tp;
if (s->parts[k].gpart != NULL) {
- s->parts[k].gpart->part = &s->parts[k];
+ s->parts[k].gpart->id_or_neg_offset = -k;
}
if (s->parts[nr_parts].gpart != NULL) {
- s->parts[nr_parts].gpart->part = &s->parts[nr_parts];
+ s->parts[nr_parts].gpart->id_or_neg_offset = -nr_parts;
}
const struct xpart txp = s->xparts[k];
s->xparts[k] = s->xparts[nr_parts];
@@ -481,11 +481,12 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
const struct gpart tp = s->gparts[k];
s->gparts[k] = s->gparts[nr_gparts];
s->gparts[nr_gparts] = tp;
- if (s->gparts[k].id > 0) {
- s->gparts[k].part->gpart = &s->gparts[k];
+ if (s->gparts[k].id_or_neg_offset <= 0) {
+ s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
- if (s->gparts[nr_gparts].id > 0) {
- s->gparts[nr_gparts].part->gpart = &s->gparts[nr_gparts];
+ if (s->gparts[nr_gparts].id_or_neg_offset <= 0) {
+ s->parts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
+ &s->gparts[nr_gparts];
}
const int t = gind[k];
gind[k] = gind[nr_gparts];
@@ -551,8 +552,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
/* Re-link the gparts. */
- for (size_t k = 0; k < nr_parts; k++)
- if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
+ part_relink_gparts(s->parts, nr_parts, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify space_sort_struct. */
@@ -600,8 +600,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
space_gparts_sort(s, gind, nr_gparts, 0, s->nr_cells - 1, verbose);
/* Re-link the parts. */
- for (int k = 0; k < nr_gparts; k++)
- if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
+ part_relink_parts(s->gparts, nr_gparts, s->parts);
/* We no longer need the indices as of here. */
free(gind);
@@ -610,21 +609,22 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Verify that the links are correct */
for (size_t k = 0; k < nr_gparts; ++k) {
- if (s->gparts[k].id > 0) {
+ if (s->gparts[k].id_or_neg_offset < 0) {
- if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+ const struct part *part = &s->parts[-s->gparts[k].id_or_neg_offset];
- if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
- s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
- s->gparts[k].x[2] != s->gparts[k].part->x[2])
+ if (part->gpart != &s->gparts[k]) error("Linking problem !");
+
+ if (s->gparts[k].x[0] != part->x[0] || s->gparts[k].x[1] != part->x[1] ||
+ s->gparts[k].x[2] != part->x[2])
error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < nr_parts; ++k) {
- if (s->parts[k].gpart != NULL) {
-
- if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
+ if (s->parts[k].gpart != NULL &&
+ s->parts[k].gpart->id_or_neg_offset != -k) {
+ error("Linking problem !");
}
}
#endif
@@ -1268,7 +1268,7 @@ void space_do_split(struct space *s, struct cell *c) {
}
/* Split the cell data. */
- cell_split(c);
+ cell_split(c, c->parts - s->parts);
/* Remove any progeny with zero parts. */
for (int k = 0; k < 8; k++)
diff --git a/src/tools.c b/src/tools.c
index 1a2b794f688047183827e5c2ed6ba80ba1339080..710efd7e4dd54784ed735492b4e1af0b466ef682 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -181,16 +181,19 @@ void pairs_single_density(double *dim, long long int pid,
}
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
+
float r2, hi, hj, hig2, hjg2, dx[3];
struct part *pi, *pj;
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = 0; j < cj->count; ++j) {
+
pj = &cj->parts[j];
/* Pairwise distance */
@@ -202,6 +205,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Hit or miss? */
if (r2 < hig2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj);
}
@@ -210,11 +214,13 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Reverse double-for loop and checks every interaction */
for (int j = 0; j < cj->count; ++j) {
+
pj = &cj->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
/* Pairwise distance */
@@ -226,6 +232,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
/* Hit or miss? */
if (r2 < hjg2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi);
}
@@ -239,11 +246,13 @@ void self_all_density(struct runner *r, struct cell *ci) {
/* Implements a double-for loop and checks every interaction */
for (int i = 0; i < ci->count; ++i) {
+
pi = &ci->parts[i];
hi = pi->h;
hig2 = hi * hi * kernel_gamma2;
for (int j = i + 1; j < ci->count; ++j) {
+
pj = &ci->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
@@ -259,12 +268,14 @@ void self_all_density(struct runner *r, struct cell *ci) {
/* Hit or miss? */
if (r2 < hig2) {
+
/* Interact */
runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj);
}
/* Hit or miss? */
if (r2 < hjg2) {
+
dxi[0] = -dxi[0];
dxi[1] = -dxi[1];
dxi[2] = -dxi[2];
@@ -277,7 +288,9 @@ void self_all_density(struct runner *r, struct cell *ci) {
}
void pairs_single_grav(double *dim, long long int pid,
- struct gpart *restrict parts, int N, int periodic) {
+ struct gpart *restrict gparts, const struct part *parts,
+ int N, int periodic) {
+
int i, k;
// int mj, mk;
// double maxratio = 1.0;
@@ -288,18 +301,20 @@ void pairs_single_grav(double *dim, long long int pid,
/* Find "our" part. */
for (k = 0; k < N; k++)
- if ((parts[k].id > 0 && parts[k].part->id == pid) || parts[k].id == -pid)
+ if ((gparts[k].id_or_neg_offset < 0 &&
+ parts[-gparts[k].id_or_neg_offset].id == pid) ||
+ gparts[k].id_or_neg_offset == pid)
break;
if (k == N) error("Part not found.");
- pi = parts[k];
+ pi = gparts[k];
pi.a_grav[0] = 0.0f;
pi.a_grav[1] = 0.0f;
pi.a_grav[2] = 0.0f;
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
- if (parts[k].id == pi.id) continue;
- pj = parts[k];
+ if (gparts[k].id_or_neg_offset == pi.id_or_neg_offset) continue;
+ pj = gparts[k];
for (i = 0; i < 3; i++) {
dx[i] = pi.x[i] - pj.x[i];
if (periodic) {
@@ -326,7 +341,8 @@ void pairs_single_grav(double *dim, long long int pid,
/* Dump the result. */
message(
"acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
- pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
+ parts[-pi.id_or_neg_offset].id, a[0], a[1], a[2], aabs[0], aabs[1],
+ aabs[2]);
}
/**
@@ -335,6 +351,7 @@ void pairs_single_grav(double *dim, long long int pid,
* @param N number of intervals in [0,1].
*/
void density_dump(int N) {
+
int k;
float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
struct part /**pi[4], *pj[4],*/ Pi[4], Pj[4];
diff --git a/src/tools.h b/src/tools.h
index 2a4462f274d8e9acd1f2d8e79996ad92c630d404..97f036994bc0d1e2b4a95d806cfdbd253664a260 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -33,7 +33,8 @@
void factor(int value, int *f1, int *f2);
void density_dump(int N);
void pairs_single_grav(double *dim, long long int pid,
- struct gpart *restrict parts, int N, int periodic);
+ struct gpart *restrict gparts, const struct part *parts,
+ int N, int periodic);
void pairs_single_density(double *dim, long long int pid,
struct part *restrict parts, int N, int periodic);