diff --git a/README b/README
index 0d658c333f1328b423851031c5b5d202f43df3c2..bec538bfc35239945d60487d66ed2e02b5d363a2 100644
--- a/README
+++ b/README
@@ -15,25 +15,26 @@ Usage: swift [OPTION]... PARAMFILE
swift_mpi [OPTION]... PARAMFILE
Valid options are:
- -a Pin runners using processor affinity
- -c Run with cosmological time integration
- -C Run with cooling
+ -a Pin runners using processor affinity.
+ -c Run with cosmological time integration.
+ -C Run with cooling.
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
Allows user to check validy of parameter and IC files as well as memory limits.
- -D Always drift all particles even the ones far from active particles.
- -e Enable floating-point exceptions (debugging mode)
- -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
+ -D Always drift all particles even the ones far from active particles. This emulates
+ Gadget-[23] and GIZMO's default behaviours.
+ -e Enable floating-point exceptions (debugging mode).
+ -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. When unset use the time_end parameter to stop.
- -s Run with SPH
- -S Run with stars
+ -s Run with hydrodynamics.
+ -S Run with stars.
-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
2: All MPI-ranks write
- -y {int} Time-step frequency at which task graphs are dumped
- -h Print this help message and exit
+ -y {int} Time-step frequency at which task graphs are dumped.
+ -h Print this help message and exit.
See the file examples/parameter_example.yml for an example of parameter file.
diff --git a/configure.ac b/configure.ac
index d5a08c7c0863bbb6dfcdbde74d68f1fe5d909907..c94f5c23b9cc5868a199ad4d177053c6f94f639d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -200,6 +200,20 @@ if test "$enable_debugging_checks" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
fi
+# Check if gravity force checks are on for some particles.
+AC_ARG_ENABLE([gravity-force-checks],
+ [AS_HELP_STRING([--enable-gravity-force-checks],
+ [Activate expensive brute-force gravity checks for a fraction 1/N of all particles @<:@N@:>@]
+ )],
+ [gravity_force_checks="$enableval"],
+ [gravity_force_checks="no"]
+)
+if test "$gravity_force_checks" == "yes"; then
+ AC_MSG_ERROR(Need to specify the fraction of particles to check when using --enable-gravity-force-checks!)
+elif test "$gravity_force_checks" != "no"; then
+ AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
+fi
+
# Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN
@@ -637,13 +651,13 @@ AC_ARG_WITH([hydro-dimension],
)
case "$with_dimension" in
1)
- AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D solver])
;;
2)
- AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D solver])
;;
3)
- AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D solver])
;;
*)
AC_MSG_ERROR([Dimensionality must be 1, 2 or 3])
@@ -824,8 +838,10 @@ AC_MSG_RESULT([
Riemann solver : $with_riemann
Cooling function : $with_cooling
External potential : $with_potential
+
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
+ Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
diff --git a/examples/main.c b/examples/main.c
index 3c1f4aab261bd7e031a42d408f4275d766948e86..fd57d632b2e649730873286f049908d3abe3178f 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -57,9 +57,9 @@ void print_help_message() {
printf(" swift_mpi [OPTION]... PARAMFILE\n\n");
printf("Valid options are:\n");
- printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
- printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
- printf(" %2s %8s %s\n", "-C", "", "Run with cooling");
+ printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity.");
+ printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration.");
+ printf(" %2s %8s %s\n", "-C", "", "Run with cooling.");
printf(
" %2s %8s %s\n", "-d", "",
"Dry run. Read the parameter file, allocate memory but does not read ");
@@ -70,29 +70,32 @@ void print_help_message() {
"Allows user to check validy of parameter and IC files as well as "
"memory limits.");
printf(" %2s %8s %s\n", "-D", "",
- "Always drift all particles even the ones far from active particles.");
+ "Always drift all particles even the ones far from active particles. "
+ "This emulates");
+ printf(" %2s %8s %s\n", "", "",
+ "Gadget-[23] and GIZMO's default behaviours.");
printf(" %2s %8s %s\n", "-e", "",
- "Enable floating-point exceptions (debugging mode)");
+ "Enable floating-point exceptions (debugging mode).");
printf(" %2s %8s %s\n", "-f", "{int}",
- "Overwrite the CPU frequency (Hz) to be used for time measurements");
+ "Overwrite the CPU frequency (Hz) to be used for time measurements.");
printf(" %2s %8s %s\n", "-g", "",
- "Run with an external gravitational potential");
- printf(" %2s %8s %s\n", "-F", "", "Run with feedback ");
- printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity");
+ "Run with an external gravitational potential.");
+ printf(" %2s %8s %s\n", "-F", "", "Run with feedback.");
+ printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity.");
printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps. When unset use the time_end "
"parameter to stop.");
- printf(" %2s %8s %s\n", "-s", "", "Run with SPH");
- printf(" %2s %8s %s\n", "-S", "", "Run with stars");
+ printf(" %2s %8s %s\n", "-s", "", "Run with hydrodynamics.");
+ printf(" %2s %8s %s\n", "-S", "", "Run with stars.");
printf(" %2s %8s %s\n", "-t", "{int}",
"The number of threads to use on each MPI rank. Defaults to 1 if not "
"specified.");
- printf(" %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity");
+ printf(" %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
printf(" %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
printf(" %2s %8s %s\n", "", "", "2: All MPI-ranks write");
printf(" %2s %8s %s\n", "-y", "{int}",
- "Time-step frequency at which task graphs are dumped");
- printf(" %2s %8s %s\n", "-h", "", "Print this help message and exit");
+ "Time-step frequency at which task graphs are dumped.");
+ printf(" %2s %8s %s\n", "-h", "", "Print this help message and exit.");
printf(
"\nSee the file parameter_example.yml for an example of "
"parameter file.\n");
@@ -302,6 +305,15 @@ int main(int argc, char *argv[]) {
message("WARNING: Debugging checks activated. Code will be slower !");
#endif
+/* Do we have gravity accuracy checks ? */
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ if (myrank == 0)
+ message(
+ "WARNING: Checking 1/%d of all gpart for gravity accuracy. Code will "
+ "be slower !",
+ SWIFT_GRAVITY_FORCE_CHECKS);
+#endif
+
/* Do we choke on FP-exceptions ? */
if (with_fp_exceptions) {
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
@@ -311,13 +323,14 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */
if (myrank == 0) {
- message("sizeof(part) is %4zi bytes.", sizeof(struct part));
- message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
- message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
- message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
- message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
- message("sizeof(task) is %4zi bytes.", sizeof(struct task));
- message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
+ message("sizeof(part) is %4zi bytes.", sizeof(struct part));
+ message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
+ message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
+ message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
+ message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
+ message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
+ message("sizeof(task) is %4zi bytes.", sizeof(struct task));
+ message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
}
/* Read the parameter file */
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index bdcd880c83eb02fdbed8d9d4da8642f3a90fc97b..7c69d19bc215472c90bb1b91f23f46702afc64f2 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -47,6 +47,7 @@ SPH:
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
max_volume_change: 2. # (Optional) Maximal allowed change of kernel volume over one time-step
+ h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
# Parameters for the self-gravity scheme
Gravity:
@@ -54,7 +55,6 @@ Gravity:
epsilon: 0.1 # Softening length (in internal units).
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
-
# Parameters related to the initial conditions
InitialConditions:
diff --git a/src/Makefile.am b/src/Makefile.am
index 2fc7ac5134c793f9476136759676a3ffaaf9816d..f2d641637554fc71af61246e434b90fa334b8b32 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -55,7 +55,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
- part_type.c xmf.c gravity_properties.c
+ part_type.c xmf.c gravity_properties.c gravity.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/cell.c b/src/cell.c
index 734b8558512bc99f60f48e7e173b57507f0ece7d..753bdd55061ebd2b2cdd2067691bd8319ca5a623 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1056,6 +1056,23 @@ void cell_check_drift_point(struct cell *c, void *data) {
#endif
}
+/**
+ * @brief Resets all the individual cell task counters to 0.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param c The #cell to reset.
+ */
+void cell_reset_task_counters(struct cell *c) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int t = 0; t < task_type_count; ++t) c->tasks_executed[t] = 0;
+ for (int t = 0; t < task_subtype_count; ++t) c->subtasks_executed[t] = 0;
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
/**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
@@ -1095,7 +1112,8 @@ int cell_are_neighbours(const struct cell *restrict ci,
*/
void cell_check_multipole(struct cell *c, void *data) {
- struct multipole ma;
+#ifdef SWIFT_DEBUG_CHECKS
+ struct gravity_tensors ma;
const double tolerance = 1e-5; /* Relative */
/* First recurse */
@@ -1106,18 +1124,22 @@ void cell_check_multipole(struct cell *c, void *data) {
if (c->gcount > 0) {
/* Brute-force calculation */
- multipole_P2M(&ma, c->gparts, c->gcount);
+ gravity_P2M(&ma, c->gparts, c->gcount);
/* Now compare the multipole expansion */
- if (!multipole_equal(&ma, c->multipole, tolerance)) {
+ if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole,
+ tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth);
message("Correct answer:");
- multipole_print(&ma);
+ gravity_multipole_print(&ma.m_pole);
message("Recursive multipole:");
- multipole_print(c->multipole);
+ gravity_multipole_print(&c->multipole->m_pole);
error("Aborting");
}
}
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
}
/**
@@ -1299,6 +1321,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+ if (c->grav_down != NULL) scheduler_activate(s, c->grav_down);
+ if (c->grav_long_range != NULL) scheduler_activate(s, c->grav_long_range);
+ if (c->grav_top_level != NULL) scheduler_activate(s, c->grav_top_level);
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
@@ -1333,6 +1358,7 @@ void cell_set_super(struct cell *c, struct cell *super) {
*/
void cell_drift_particles(struct cell *c, const struct engine *e) {
+ const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
const integertime_t ti_old = c->ti_old;
const integertime_t ti_current = e->ti_current;
@@ -1343,7 +1369,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
- float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
+ float dx_max = 0.f, dx2_max = 0.f, cell_h_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
@@ -1357,7 +1383,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
cell_drift_particles(cp, e);
dx_max = max(dx_max, cp->dx_max);
- h_max = max(h_max, cp->h_max);
+ cell_h_max = max(cell_h_max, cp->h_max);
}
} else if (ti_current > ti_old) {
@@ -1376,7 +1402,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
}
/* Loop over all the gas particles in the cell */
@@ -1390,14 +1416,17 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
+ /* Limit h to within the allowed range */
+ p->h = min(p->h, hydro_h_max);
+
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
/* Maximal smoothing length */
- h_max = (h_max > p->h) ? h_max : p->h;
+ cell_h_max = max(cell_h_max, p->h);
}
/* Loop over all the star particles in the cell */
@@ -1418,12 +1447,12 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
} else {
- h_max = c->h_max;
+ cell_h_max = c->h_max;
dx_max = c->dx_max;
}
/* Store the values */
- c->h_max = h_max;
+ c->h_max = cell_h_max;
c->dx_max = dx_max;
/* Update the time of the last drift */
@@ -1436,7 +1465,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
-void cell_drift_multipole(struct cell *c, const struct engine *e) {
+void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
const double timeBase = e->timeBase;
const integertime_t ti_old_multipole = c->ti_old_multipole;
@@ -1458,13 +1487,26 @@ void cell_drift_multipole(struct cell *c, const struct engine *e) {
} else if (ti_current > ti_old_multipole) {
/* Drift the multipole */
- multipole_drift(c->multipole, dt);
+ gravity_multipole_drift(c->multipole, dt);
}
/* Update the time of the last drift */
c->ti_old_multipole = ti_current;
}
+/**
+ * @brief Drifts the multipole of a cell to the current time.
+ *
+ * Only drifts the multipole at this level. Multipoles deeper in the
+ * tree are not updated.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_multipole(struct cell *c, const struct engine *e) {
+ error("To be implemented");
+}
+
/**
* @brief Recursively checks that all particles in a cell have a time-step
*/
diff --git a/src/cell.h b/src/cell.h
index dc0e8886e0ae202df555049925036804b37b6f6a..b9fbb6abffe427b81daec033447fbf35804c242d 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -104,7 +104,7 @@ struct cell {
double h_max;
/*! This cell's multipole. */
- struct multipole *multipole;
+ struct gravity_tensors *multipole;
/*! Linking pointer for "memory management". */
struct cell *next;
@@ -170,7 +170,10 @@ struct cell {
struct task *timestep;
/*! Task constructing the multipole from the particles */
- struct task *grav_up;
+ struct task *grav_top_level;
+
+ /*! Task constructing the multipole from the particles */
+ struct task *grav_long_range;
/*! Task propagating the multipole to the particles */
struct task *grav_down;
@@ -308,6 +311,14 @@ struct cell {
/*! The maximal depth of this cell and its progenies */
char maxdepth;
+#ifdef SWIFT_DEBUG_CHECKS
+ /*! The list of tasks that have been executed on this cell */
+ char tasks_executed[64];
+
+ /*! The list of sub-tasks that have been executed on this cell */
+ char subtasks_executed[64];
+#endif
+
} SWIFT_STRUCT_ALIGN;
/* Convert cell location to ID. */
@@ -342,11 +353,13 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
void cell_check_drift_point(struct cell *c, void *data);
+void cell_reset_task_counters(struct cell *c);
int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_particles(struct cell *c, const struct engine *e);
void cell_drift_multipole(struct cell *c, const struct engine *e);
+void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
#endif /* SWIFT_CELL_H */
diff --git a/src/const.h b/src/const.h
index 010cad4167c28c93a69d9e23f7dc8612462c20bf..88c1a1af1cfc36cc401fdfea0b077f79fcd13bc0 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,9 +41,6 @@
/* Self gravity stuff. */
#define const_gravity_multipole_order 1
-#define const_gravity_a_smooth 1.25f
-#define const_gravity_r_cut 4.5f
-#define const_gravity_eta 0.025f
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
diff --git a/src/engine.c b/src/engine.c
index 425e671d97fbd2e47c4f1d781737ecf943c84ac0..facffa389e98fdeeebda222f04b2953c00b97373 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -55,6 +55,7 @@
#include "cycle.h"
#include "debug.h"
#include "error.h"
+#include "gravity.h"
#include "hydro.h"
#include "minmax.h"
#include "parallel_io.h"
@@ -132,6 +133,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
const int is_hydro = (e->policy & engine_policy_hydro);
+ const int is_self_gravity = (e->policy & engine_policy_self_gravity);
const int is_with_cooling = (e->policy & engine_policy_cooling);
const int is_with_sourceterms = (e->policy & engine_policy_sourceterms);
@@ -165,6 +167,27 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->drift, c->init);
+ if (is_self_gravity) {
+
+ /* Gravity non-neighbouring pm calculations */
+ c->grav_long_range = scheduler_addtask(
+ s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL, 0);
+
+ /* Gravity top-level periodic calculation */
+ c->grav_top_level = scheduler_addtask(
+ s, task_type_grav_top_level, task_subtype_none, 0, 0, c, NULL, 0);
+
+ /* Gravity recursive down-pass */
+ c->grav_down = scheduler_addtask(s, task_type_grav_down,
+ task_subtype_none, 0, 0, c, NULL, 0);
+
+ scheduler_addunlock(s, c->init, c->grav_long_range);
+ scheduler_addunlock(s, c->init, c->grav_top_level);
+ scheduler_addunlock(s, c->grav_long_range, c->grav_down);
+ scheduler_addunlock(s, c->grav_top_level, c->grav_down);
+ scheduler_addunlock(s, c->grav_down, c->kick2);
+ }
+
/* Generate the ghost task. */
if (is_hydro)
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
@@ -871,7 +894,7 @@ void engine_addtasks_grav(struct engine *e, struct cell *c, struct task *up,
struct task *down) {
/* Link the tasks to this cell. */
- c->grav_up = up;
+ // c->grav_up = up;
c->grav_down = down;
/* Recurse? */
@@ -1550,7 +1573,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
*
* @param e The #engine.
*/
-void engine_make_gravity_tasks(struct engine *e) {
+void engine_make_self_gravity_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
@@ -1572,10 +1595,6 @@ void engine_make_gravity_tasks(struct engine *e) {
scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL,
0);
- /* Let's also build a task for all the non-neighbouring pm calculations */
- scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, 0, 0, ci,
- NULL, 0);
-
for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
struct cell *cj = &cells[cjd];
@@ -1777,25 +1796,20 @@ void engine_count_and_link_tasks(struct engine *e) {
}
}
-/*
+/**
* @brief Creates the dependency network for the gravity tasks of a given cell.
*
* @param sched The #scheduler.
* @param gravity The gravity task to link.
* @param c The cell.
*/
-/* static inline void engine_make_gravity_dependencies(struct scheduler *sched,
- */
-/* struct task *gravity, */
-/* struct cell *c) { */
-
-/* /\* init --> gravity --> kick *\/ */
-/* scheduler_addunlock(sched, c->super->init, gravity); */
-/* scheduler_addunlock(sched, gravity, c->super->kick2); */
+static inline void engine_make_self_gravity_dependencies(
+ struct scheduler *sched, struct task *gravity, struct cell *c) {
-/* /\* grav_up --> gravity ( --> kick) *\/ */
-/* scheduler_addunlock(sched, c->super->grav_up, gravity); */
-/* } */
+ /* init --> gravity --> grav_down --> kick */
+ scheduler_addunlock(sched, c->super->init, gravity);
+ scheduler_addunlock(sched, gravity, c->super->grav_down);
+}
/**
* @brief Creates the dependency network for the external gravity tasks of a
@@ -1829,12 +1843,11 @@ void engine_link_gravity_tasks(struct engine *e) {
/* Get a pointer to the task. */
struct task *t = &sched->tasks[k];
- /* /\* Self-interaction for self-gravity? *\/ */
- /* if (t->type == task_type_self && t->subtype == task_subtype_grav) { */
-
- /* engine_make_gravity_dependencies(sched, t, t->ci); */
+ /* Self-interaction for self-gravity? */
+ if (t->type == task_type_self && t->subtype == task_subtype_grav) {
- /* } */
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
+ }
/* Self-interaction for external gravity ? */
if (t->type == task_type_self && t->subtype == task_subtype_external_grav) {
@@ -1843,31 +1856,28 @@ void engine_link_gravity_tasks(struct engine *e) {
}
- /* /\* Otherwise, pair interaction? *\/ */
- /* else if (t->type == task_type_pair && t->subtype == task_subtype_grav)
- * {
- */
+ /* Otherwise, pair interaction? */
+ else if (t->type == task_type_pair && t->subtype == task_subtype_grav) {
- /* if (t->ci->nodeID == nodeID) { */
+ if (t->ci->nodeID == nodeID) {
- /* engine_make_gravity_dependencies(sched, t, t->ci); */
- /* } */
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
+ }
- /* if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) { */
+ if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
- /* engine_make_gravity_dependencies(sched, t, t->cj); */
- /* } */
+ engine_make_self_gravity_dependencies(sched, t, t->cj);
+ }
- /* } */
+ }
- /* /\* Otherwise, sub-self interaction? *\/ */
- /* else if (t->type == task_type_sub_self && t->subtype ==
- * task_subtype_grav) { */
+ /* Otherwise, sub-self interaction? */
+ else if (t->type == task_type_sub_self && t->subtype == task_subtype_grav) {
- /* if (t->ci->nodeID == nodeID) { */
- /* engine_make_gravity_dependencies(sched, t, t->ci); */
- /* } */
- /* } */
+ if (t->ci->nodeID == nodeID) {
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
+ }
+ }
/* Sub-self-interaction for external gravity ? */
else if (t->type == task_type_sub_self &&
@@ -1878,19 +1888,18 @@ void engine_link_gravity_tasks(struct engine *e) {
}
}
- /* /\* Otherwise, sub-pair interaction? *\/ */
- /* else if (t->type == task_type_sub_pair && t->subtype ==
- * task_subtype_grav) { */
+ /* Otherwise, sub-pair interaction? */
+ else if (t->type == task_type_sub_pair && t->subtype == task_subtype_grav) {
- /* if (t->ci->nodeID == nodeID) { */
+ if (t->ci->nodeID == nodeID) {
- /* engine_make_gravity_dependencies(sched, t, t->ci); */
- /* } */
- /* if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) { */
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
+ }
+ if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
- /* engine_make_gravity_dependencies(sched, t, t->cj); */
- /* } */
- /* } */
+ engine_make_self_gravity_dependencies(sched, t, t->cj);
+ }
+ }
}
}
@@ -2180,17 +2189,11 @@ void engine_make_gravityrecursive_tasks(struct engine *e) {
/* if (cells[k].nodeID == nodeID && cells[k].gcount > 0) { */
/* /\* Create tasks at top level. *\/ */
- /* struct task *up = */
- /* scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0,
- * 0, */
- /* &cells[k], NULL, 0); */
-
+ /* struct task *up = NULL; */
/* struct task *down = NULL; */
- /* /\* struct task *down = *\/ */
- /* /\* scheduler_addtask(sched, task_type_grav_down,
- * task_subtype_none, 0, */
- /* * 0, *\/ */
- /* /\* &cells[k], NULL, 0); *\/ */
+ /* /\* scheduler_addtask(sched, task_type_grav_down,
+ * task_subtype_none, 0, 0, *\/ */
+ /* /\* &cells[k], NULL, 0); *\/ */
/* /\* Push tasks down the cell hierarchy. *\/ */
/* engine_addtasks_grav(e, &cells[k], up, down); */
@@ -2217,8 +2220,8 @@ void engine_maketasks(struct engine *e) {
/* Construct the firt hydro loop over neighbours */
if (e->policy & engine_policy_hydro) engine_make_hydroloop_tasks(e);
- /* Add the gravity mm tasks. */
- if (e->policy & engine_policy_self_gravity) engine_make_gravity_tasks(e);
+ /* Add the self gravity tasks. */
+ if (e->policy & engine_policy_self_gravity) engine_make_self_gravity_tasks(e);
/* Add the external gravity tasks. */
if (e->policy & engine_policy_external_gravity)
@@ -2461,6 +2464,13 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
+ /* Gravity ? */
+ else if (t->type == task_type_grav_down ||
+ t->type == task_type_grav_long_range ||
+ t->type == task_type_grav_top_level) {
+ if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
+ }
+
/* Time-step? */
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
@@ -2531,6 +2541,7 @@ void engine_print_task_counts(struct engine *e) {
fflush(stdout);
message("nr_parts = %zu.", e->s->nr_parts);
message("nr_gparts = %zu.", e->s->nr_gparts);
+ message("nr_sparts = %zu.", e->s->nr_sparts);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2839,8 +2850,10 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift || t->type == task_type_kick1 ||
t->type == task_type_kick2 || t->type == task_type_timestep ||
- t->subtype == task_subtype_force || t->type == task_type_cooling ||
- t->type == task_type_sourceterms)
+ t->subtype == task_subtype_force || t->subtype == task_subtype_grav ||
+ t->type == task_type_grav_long_range ||
+ t->type == task_type_grav_top_level || t->type == task_type_grav_down ||
+ t->type == task_type_cooling || t->type == task_type_sourceterms)
t->skip = 1;
}
}
@@ -2874,6 +2887,11 @@ void engine_launch(struct engine *e, int nr_runners) {
const ticks tic = getticks();
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Re-set all the cell task counters to 0 */
+ space_reset_task_counters(e->s);
+#endif
+
/* Prepare the scheduler. */
atomic_inc(&e->sched.waiting);
@@ -2952,6 +2970,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
}
}
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Let's store the total mass in the system for future checks */
+ e->s->total_mass = 0.;
+ for (size_t i = 0; i < s->nr_gparts; ++i)
+ e->s->total_mass += s->gparts[i].mass;
+#ifdef WITH_MPI
+ if (MPI_Allreduce(MPI_IN_PLACE, &e->s->total_mass, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to all-reduce total mass in the system.");
+#endif
+ message("Total mass in the system: %e", e->s->total_mass);
+#endif
+
/* Now time to get ready for the first time-step */
if (e->nodeID == 0) message("Running initial fake time-step.");
@@ -3030,14 +3061,27 @@ void engine_step(struct engine *e) {
if (e->step % 100 == 2) e->forcerepart = 1;
#endif
+ /* Are we drifting everything (a la Gadget/GIZMO) ? */
+ if (e->policy & engine_policy_drift_all) engine_drift_all(e);
+
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Run the brute-force gravity calculation for some gparts */
+ gravity_exact_force_compute(e->s, e);
+#endif
+
/* Start all the tasks. */
TIMER_TIC;
engine_launch(e, e->nr_threads);
TIMER_TOC(timer_runners);
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Check the accuracy of the gravity calculation */
+ gravity_exact_force_check(e->s, e, 1e-1);
+#endif
+
/* Collect the values of rebuild from all nodes. */
#ifdef WITH_MPI
int buff = 0;
@@ -3136,7 +3180,8 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
cell_drift_particles(c, e);
/* Drift the multipole */
- if (e->policy & engine_policy_self_gravity) cell_drift_multipole(c, e);
+ if (e->policy & engine_policy_self_gravity)
+ cell_drift_all_multipoles(c, e);
}
}
}
@@ -3149,6 +3194,11 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
void engine_drift_all(struct engine *e) {
const ticks tic = getticks();
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (e->nodeID == 0) message("Drifting all");
+#endif
+
threadpool_map(&e->threadpool, engine_do_drift_all_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
diff --git a/src/gravity.c b/src/gravity.c
new file mode 100644
index 0000000000000000000000000000000000000000..86f9fa82e3eb693cc3c051420fb9c7bff277eb9f
--- /dev/null
+++ b/src/gravity.c
@@ -0,0 +1,187 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "gravity.h"
+
+/* Local headers. */
+#include "active.h"
+#include "error.h"
+
+/**
+ * @brief Run a brute-force gravity calculation for a subset of particles.
+ *
+ * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will get their forces
+ * computed.
+ *
+ * @param s The #space to use.
+ * @param e The #engine (to access the current time).
+ */
+void gravity_exact_force_compute(struct space *s, const struct engine *e) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ const ticks tic = getticks();
+ const double const_G = e->physical_constants->const_newton_G;
+ int counter = 0;
+
+ for (size_t i = 0; i < s->nr_gparts; ++i) {
+
+ struct gpart *gpi = &s->gparts[i];
+
+ /* Is the particle active and part of the subset to be tested ? */
+ if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
+ gpart_is_active(gpi, e)) {
+
+ /* Be ready for the calculation */
+ gpi->a_grav[0] = 0.f;
+ gpi->a_grav[1] = 0.f;
+ gpi->a_grav[2] = 0.f;
+
+ /* Interact it with all other particles in the space.*/
+ for (size_t j = 0; j < s->nr_gparts; ++j) {
+
+ /* No self interaction */
+ if (i == j) continue;
+
+ struct gpart *gpj = &s->gparts[j];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {gpi->x[0] - gpj->x[0], // x
+ gpi->x[1] - gpj->x[1], // y
+ gpi->x[2] - gpj->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ runner_iact_grav_pp_nonsym(0.f, r2, dx, gpi, gpj);
+ }
+
+ /* Finish the calculation */
+ gravity_end_force(gpi, const_G);
+
+ /* Store the exact answer */
+ gpi->a_grav_exact[0] = gpi->a_grav[0];
+ gpi->a_grav_exact[1] = gpi->a_grav[1];
+ gpi->a_grav_exact[2] = gpi->a_grav[2];
+
+ /* Restore everything */
+ gpi->a_grav[0] = 0.f;
+ gpi->a_grav[1] = 0.f;
+ gpi->a_grav[2] = 0.f;
+
+ counter++;
+ }
+ }
+
+ message("Computed exact gravity for %d gparts.", counter);
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+#else
+ error("Gravity checking function called without the corresponding flag.");
+#endif
+}
+
+/**
+ * @brief Check the accuracy of the gravity calculation by comparing the
+ * accelerations
+ * to the brute-force computed ones.
+ *
+ * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will be checked.
+ *
+ * @param s The #space to use.
+ * @param e The #engine (to access the current time).
+ * @param rel_tol The maximal relative error. Will call error() if one particle
+ * has a larger error.
+ */
+void gravity_exact_force_check(struct space *s, const struct engine *e,
+ float rel_tol) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ const double const_G = e->physical_constants->const_newton_G;
+
+ int counter = 0;
+
+ /* Some accumulators */
+ float err_rel[3];
+ float err_rel_max[3] = {0.f, 0.f, 0.f};
+ float err_rel_min[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
+ float err_rel_mean[3] = {0.f, 0.f, 0.f};
+ float err_rel_mean2[3] = {0.f, 0.f, 0.f};
+ float err_rel_std[3] = {0.f, 0.f, 0.f};
+
+ for (size_t i = 0; i < s->nr_gparts; ++i) {
+
+ struct gpart *gpi = &s->gparts[i];
+
+ /* Is the particle was active and part of the subset to be tested ? */
+ if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
+ gpart_is_starting(gpi, e)) {
+
+ /* Compute relative error */
+ for (int k = 0; k < 3; ++k)
+ if (fabsf(gpi->a_grav_exact[k]) > FLT_EPSILON * const_G)
+ err_rel[k] = (gpi->a_grav[k] - gpi->a_grav_exact[k]) /
+ fabsf(gpi->a_grav_exact[k]);
+ else
+ err_rel[k] = 0.f;
+
+ /* Check that we are not above tolerance */
+ if (fabsf(err_rel[0]) > rel_tol || fabsf(err_rel[1]) > rel_tol ||
+ fabsf(err_rel[2]) > rel_tol)
+ error("Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e]",
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+
+ /* Construct some statistics */
+ for (int k = 0; k < 3; ++k) {
+ err_rel_max[k] = max(err_rel_max[k], fabsf(err_rel[k]));
+ err_rel_min[k] = min(err_rel_min[k], fabsf(err_rel[k]));
+ err_rel_mean[k] += err_rel[k];
+ err_rel_mean2[k] += err_rel[k] * err_rel[k];
+ }
+
+ counter++;
+ }
+ }
+
+ /* Final operation on the stats */
+ if (counter > 0) {
+ for (int k = 0; k < 3; ++k) {
+ err_rel_mean[k] /= counter;
+ err_rel_mean2[k] /= counter;
+ err_rel_std[k] =
+ sqrtf(err_rel_mean2[k] - err_rel_mean[k] * err_rel_mean[k]);
+ }
+ }
+
+ /* Report on the findings */
+ message("Checked gravity for %d gparts.", counter);
+ for (int k = 0; k < 3; ++k)
+ message("Error on a_grav[%d]: min=%e max=%e mean=%e std=%e", k,
+ err_rel_min[k], err_rel_max[k], err_rel_mean[k], err_rel_std[k]);
+
+#else
+ error("Gravity checking function called without the corresponding flag.");
+#endif
+}
diff --git a/src/gravity.h b/src/gravity.h
index 6edcd59d9955029217364d4fc67874216b0e24dc..00b930c00fb2558f274feb2991b78e96dc8b990b 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -22,9 +22,20 @@
/* Config parameters. */
#include "../config.h"
+/* Local headers. */
+#include "const.h"
+#include "engine.h"
+#include "inline.h"
+#include "part.h"
+#include "space.h"
+
/* So far only one model here */
/* Straight-forward import */
#include "./gravity/Default/gravity.h"
#include "./gravity/Default/gravity_iact.h"
+void gravity_exact_force_compute(struct space *s, const struct engine *e);
+void gravity_exact_force_check(struct space *s, const struct engine *e,
+ float rel_tol);
+
#endif
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 3f97070fbb80c8920a929e6df8d93def6ac6041a..93a65e2f5a70e09ad14280cf9c334753359fb8b5 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -21,15 +21,18 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
+#include "gravity_properties.h"
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param gp Pointer to the g-particle data.
+ * @param grav_props Constants used in the gravity scheme.
*/
__attribute__((always_inline)) INLINE static float
-gravity_compute_timestep_self(const struct gpart* const gp) {
+gravity_compute_timestep_self(const struct gpart* const gp,
+ const struct gravity_props* restrict grav_props) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
@@ -37,7 +40,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
- const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
+ const float dt = sqrtf(2.f * grav_props->eta * gp->epsilon / ac);
return dt;
}
@@ -57,6 +60,10 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
gp->a_grav[0] = 0.f;
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->mass_interacted = 0.;
+#endif
}
/**
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index f484b13663059fa5f4f822aa78748fe4ef9d5926..00ae0f5b05cd95750c34b60e2353a9fc1d0a5c32 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -52,6 +52,9 @@ struct gpart {
#ifdef SWIFT_DEBUG_CHECKS
+ /* Total mass this gpart interacted with */
+ double mass_interacted;
+
/* Time of the last drift */
integertime_t ti_drift;
@@ -60,6 +63,13 @@ struct gpart {
#endif
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ /* Brute-force particle acceleration. */
+ float a_grav_exact[3];
+
+#endif
+
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 8db23a3de0123400fe691b0d60f076929e1b0b48..46785b4b2d5b958f6db3bd9813d139575217d6fe 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -34,6 +34,7 @@
#define hydro_props_default_max_iterations 30
#define hydro_props_default_volume_change 2.0f
+#define hydro_props_default_h_max FLT_MAX
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
@@ -43,7 +44,11 @@ void hydro_props_init(struct hydro_props *p,
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
- /* Ghost stuff */
+ /* Maximal smoothing length */
+ p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
+ hydro_props_default_h_max);
+
+ /* Number of iterations to converge h */
p->max_smoothing_iterations = parser_get_opt_param_int(
params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations);
@@ -81,6 +86,9 @@ void hydro_props_print(const struct hydro_props *p) {
"(max|dlog(h)/dt|=%f).",
pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
+ if (p->h_max != hydro_props_default_h_max)
+ message("Maximal smoothing length allowed: %.4f", p->h_max);
+
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
p->max_smoothing_iterations, hydro_props_default_max_iterations);
@@ -96,6 +104,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
+ io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
p->log_max_h_change);
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 859ae10f18a71ba920f44b8c59e63abd57cd7c92..716c4c060c21eb95d05f9d50e13d4681a958a6fd 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -40,7 +40,10 @@ struct hydro_props {
float target_neighbours;
float delta_neighbours;
- /* Kernel properties */
+ /* Maximal smoothing length */
+ float h_max;
+
+ /* Number of iterations to converge h */
int max_smoothing_iterations;
/* Time integration properties */
diff --git a/src/multipole.h b/src/multipole.h
index 9521bd51db49fd3721bd1c6ddf1c5eda1dd4c0ce..95b7450cd9c7aa60dba5c43ebd3c77e6316845a5 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -25,6 +25,7 @@
#include <string.h>
/* Includes. */
+//#include "active.h"
#include "align.h"
#include "const.h"
#include "error.h"
@@ -36,65 +37,60 @@
#define multipole_align 128
-/**
- * @brief The multipole expansion of a mass distribution.
- */
-struct multipole {
+struct acc_tensor {
- union {
+ double F_000;
+};
- /*! Linking pointer for "memory management". */
- struct multipole *next;
+struct pot_tensor {
- /*! The actual content */
- struct {
+ double F_000;
+};
- /*! Multipole mass */
- float mass;
+struct multipole {
- /*! Centre of mass of the matter dsitribution */
- double CoM[3];
+ float mass;
- /*! Bulk velocity */
- float vel[3];
- };
- };
-} SWIFT_STRUCT_ALIGN;
+ /*! Bulk velocity */
+ float vel[3];
+};
-struct acc_tensor {
+/**
+ * @brief The multipole expansion of a mass distribution.
+ */
+struct gravity_tensors {
- double F_000;
-};
+ /*! Linking pointer for "memory management". */
+ struct gravity_tensors *next;
-struct pot_tensor {
+ /*! Centre of mass of the matter dsitribution */
+ double CoM[3];
- double F_000;
-};
+ /*! The actual content */
+ struct {
-struct field_tensors {
+ /*! Multipole mass */
+ struct multipole m_pole;
- union {
+ /*! Field tensor for acceleration along x */
+ struct acc_tensor a_x;
- /*! Linking pointer for "memory management". */
- struct field_tensors *next;
+ /*! Field tensor for acceleration along y */
+ struct acc_tensor a_y;
- /*! The actual content */
- struct {
+ /*! Field tensor for acceleration along z */
+ struct acc_tensor a_z;
- /*! Field tensor for acceleration along x */
- struct acc_tensor a_x;
+ /*! Field tensor for the potential */
+ struct pot_tensor pot;
- /*! Field tensor for acceleration along y */
- struct acc_tensor a_y;
+#ifdef SWIFT_DEBUG_CHECKS
- /*! Field tensor for acceleration along z */
- struct acc_tensor a_z;
+ /* Total mass this gpart interacted with */
+ double mass_interacted;
- /*! Field tensor for the potential */
- struct pot_tensor pot;
- };
+#endif
};
-
} SWIFT_STRUCT_ALIGN;
/**
@@ -102,10 +98,22 @@ struct field_tensors {
*
* @param m The #multipole.
*/
-INLINE static void multipole_init(struct multipole *m) {
+INLINE static void gravity_reset(struct gravity_tensors *m) {
/* Just bzero the struct. */
- bzero(m, sizeof(struct multipole));
+ bzero(m, sizeof(struct gravity_tensors));
+}
+
+INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
+
+ bzero(&m->a_x, sizeof(struct acc_tensor));
+ bzero(&m->a_y, sizeof(struct acc_tensor));
+ bzero(&m->a_z, sizeof(struct acc_tensor));
+ bzero(&m->pot, sizeof(struct pot_tensor));
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m->mass_interacted = 0.;
+#endif
}
/**
@@ -115,9 +123,9 @@ INLINE static void multipole_init(struct multipole *m) {
*
* @param m The #multipole to print.
*/
-INLINE static void multipole_print(const struct multipole *m) {
+INLINE static void gravity_multipole_print(const struct multipole *m) {
- printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
+ // printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
printf("Mass= %12.5e\n", m->mass);
printf("Vel= [%12.5e %12.5e %12.5e\n", m->vel[0], m->vel[1], m->vel[2]);
}
@@ -128,8 +136,8 @@ INLINE static void multipole_print(const struct multipole *m) {
* @param ma The multipole to add to.
* @param mb The multipole to add.
*/
-INLINE static void multipole_add(struct multipole *ma,
- const struct multipole *mb) {
+INLINE static void gravity_multipole_add(struct multipole *ma,
+ const struct multipole *mb) {
const float mass = ma->mass + mb->mass;
const float imass = 1.f / mass;
@@ -151,17 +159,20 @@ INLINE static void multipole_add(struct multipole *ma,
* @param tolerance The maximal allowed difference for the fields.
* @return 1 if the multipoles are equal 0 otherwise.
*/
-INLINE static int multipole_equal(const struct multipole *ma,
- const struct multipole *mb,
- double tolerance) {
+INLINE static int gravity_multipole_equal(const struct multipole *ma,
+ const struct multipole *mb,
+ double tolerance) {
/* Check CoM */
- if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) > tolerance)
- return 0;
- if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) > tolerance)
- return 0;
- if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) > tolerance)
- return 0;
+ /* if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) >
+ * tolerance) */
+ /* return 0; */
+ /* if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) >
+ * tolerance) */
+ /* return 0; */
+ /* if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) >
+ * tolerance) */
+ /* return 0; */
/* Check bulk velocity (if non-zero)*/
if (fabsf(ma->vel[0] + mb->vel[0]) > 0.f &&
@@ -191,12 +202,13 @@ INLINE static int multipole_equal(const struct multipole *ma,
* @param m The #multipole.
* @param dt The drift time-step.
*/
-INLINE static void multipole_drift(struct multipole *m, double dt) {
+INLINE static void gravity_multipole_drift(struct gravity_tensors *m,
+ double dt) {
/* Move the whole thing according to bulk motion */
- m->CoM[0] += m->vel[0];
- m->CoM[1] += m->vel[1];
- m->CoM[2] += m->vel[2];
+ m->CoM[0] += m->m_pole.vel[0];
+ m->CoM[1] += m->m_pole.vel[1];
+ m->CoM[2] += m->m_pole.vel[2];
}
/**
@@ -207,8 +219,8 @@ INLINE static void multipole_drift(struct multipole *m, double dt) {
* @param gparts_j The #gpart that source the gravity field.
* @param gcount_j The number of sources.
*/
-INLINE static void multipole_P2P(struct gpart *gparts_i, int gcount_i,
- const struct gpart *gparts_j, int gcount_j) {}
+INLINE static void gravity_P2P(struct gpart *gparts_i, int gcount_i,
+ const struct gpart *gparts_j, int gcount_j) {}
/**
* @brief Constructs the #multipole of a bunch of particles around their
@@ -220,8 +232,8 @@ INLINE static void multipole_P2P(struct gpart *gparts_i, int gcount_i,
* @param gparts The #gpart.
* @param gcount The number of particles.
*/
-INLINE static void multipole_P2M(struct multipole *m,
- const struct gpart *gparts, int gcount) {
+INLINE static void gravity_P2M(struct gravity_tensors *m,
+ const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
#error "Implementation of P2M kernel missing for this order."
@@ -248,13 +260,13 @@ INLINE static void multipole_P2M(struct multipole *m,
const double imass = 1.0 / mass;
/* Store the data on the multipole. */
- m->mass = mass;
+ m->m_pole.mass = mass;
m->CoM[0] = com[0] * imass;
m->CoM[1] = com[1] * imass;
m->CoM[2] = com[2] * imass;
- m->vel[0] = vel[0] * imass;
- m->vel[1] = vel[1] * imass;
- m->vel[2] = vel[2] * imass;
+ m->m_pole.vel[0] = vel[0] * imass;
+ m->m_pole.vel[1] = vel[1] * imass;
+ m->m_pole.vel[2] = vel[2] * imass;
}
/**
@@ -268,10 +280,10 @@ INLINE static void multipole_P2M(struct multipole *m,
* @param pos_b The current postion of the multipole to shift.
* @param periodic Is the calculation periodic ?
*/
-INLINE static void multipole_M2M(struct multipole *m_a,
- const struct multipole *m_b,
- const double pos_a[3], const double pos_b[3],
- int periodic) {
+INLINE static void gravity_M2M(struct multipole *m_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
m_a->mass = m_b->mass;
@@ -290,38 +302,79 @@ INLINE static void multipole_M2M(struct multipole *m_a,
* @param pos_b The position of the multipole.
* @param periodic Is the calculation periodic ?
*/
-INLINE static void multipole_M2L(struct field_tensors *l_a,
- const struct multipole m_b,
- const double pos_a[3], const double pos_b[3],
- int periodic) {
-
- /* double dx, dy, dz; */
- /* if (periodic) { */
- /* dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.); */
- /* dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.); */
- /* dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.); */
- /* } else { */
- /* dx = pos_a[0] - pos_b[0]; */
- /* dy = pos_a[1] - pos_b[1]; */
- /* dz = pos_a[2] - pos_b[2]; */
- /* } */
- /* const double r2 = dx * dx + dy * dy + dz * dz; */
-
- /* const double r_inv = 1. / sqrt(r2); */
-
- /* /\* 1st order multipole term *\/ */
- /* l_a->x.F_000 = D_100(dx, dy, dz, r_inv); */
- /* l_a->y.F_000 = D_010(dx, dy, dz, r_inv); */
- /* l_a->z.F_000 = D_001(dx, dy, dz, r_inv); */
+INLINE static void gravity_M2L(struct gravity_tensors *l_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
+
+ double dx, dy, dz;
+ if (periodic) {
+ dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.);
+ dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.);
+ dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.);
+ } else {
+ dx = pos_a[0] - pos_b[0];
+ dy = pos_a[1] - pos_b[1];
+ dz = pos_a[2] - pos_b[2];
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ const double r_inv = 1. / sqrt(r2);
+
+ /* 1st order multipole term */
+ l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ l_a->mass_interacted += m_b->mass;
+#endif
+}
+
+/**
+ * @brief Creates a copy of #acc_tensor shifted to a new location.
+ *
+ * Corresponds to equation (28e).
+ *
+ * @param l_a The #acc_tensor copy (content will be overwritten).
+ * @param l_b The #acc_tensor to shift.
+ * @param pos_a The position to which m_b will be shifted.
+ * @param pos_b The current postion of the multipole to shift.
+ * @param periodic Is the calculation periodic ?
+ */
+INLINE static void gravity_L2L(struct gravity_tensors *l_a,
+ const struct gravity_tensors *l_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {}
+
+/**
+ * @brief Applies the #acc_tensor to a set of #gpart.
+ *
+ * Corresponds to equation (28a).
+ */
+INLINE static void gravity_L2P(const struct gravity_tensors *l,
+ struct gpart *gparts, int gcount) {
+
+ for (int i = 0; i < gcount; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ struct gpart *gp = &gparts[i];
+
+// if(gpart_is_active(gp, e)){
+
+ gp->mass_interacted += l->mass_interacted;
+#endif
+ //}
+ }
}
#if 0
/* Multipole function prototypes. */
-void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
-void multipole_init(struct multipole *m, const struct gpart *gparts,
+void multipole_add(struct gravity_tensors *m_sum, const struct gravity_tensors *m_term);
+void multipole_init(struct gravity_tensors *m, const struct gpart *gparts,
int gcount);
-void multipole_reset(struct multipole *m);
+void multipole_reset(struct gravity_tensors *m);
/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
/* double *shift); */
@@ -336,7 +389,7 @@ void multipole_reset(struct multipole *m);
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
- struct multipole *ma, struct multipole *mb, double *shift) {
+ struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
@@ -378,7 +431,7 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mm(
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
- struct multipole *m, struct gpart *p, double *shift) {
+ struct gravity_tensors *m, struct gpart *p, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
diff --git a/src/part.c b/src/part.c
index 98033097d749d1aff927e0781020c2ea25d74ffd..712a13d42a629355290bcf89f31b94fb9670a215 100644
--- a/src/part.c
+++ b/src/part.c
@@ -267,8 +267,9 @@ void part_create_mpi_types() {
MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for sparts.");
}
- if (MPI_Type_contiguous(sizeof(struct multipole) / sizeof(unsigned char),
- MPI_BYTE, &multipole_mpi_type) != MPI_SUCCESS ||
+ if (MPI_Type_contiguous(
+ sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
+ &multipole_mpi_type) != MPI_SUCCESS ||
MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for multipole.");
}
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index 160b74cae82f8dba584e15eac743e6ddc2002c0a..adea9d912056fd07e134fb98a7603030e897ec7a 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -24,6 +24,7 @@
#include "../config.h"
/* Some standard headers. */
+#include <float.h>
#include <math.h>
/* Local includes. */
@@ -84,7 +85,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
- return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ if (fabsf(dota_2) > 0.f)
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ else
+ return FLT_MAX;
}
/**
diff --git a/src/runner.c b/src/runner.c
index f0970ae701ad0bd888cefc20b26d73441e21d298..521914bfdecb18f6c5e4e0a55ec2f0e759424797 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -113,7 +113,7 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* Import the gravity loop functions. */
#include "runner_doiact_fft.h"
-//#include "runner_doiact_grav.h"
+#include "runner_doiact_grav.h"
/**
* @brief Perform source terms
@@ -495,6 +495,10 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
+ /* Reset the gravity acceleration tensors */
+ if (e->policy & engine_policy_self_gravity)
+ gravity_field_tensor_init(c->multipole);
+
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
@@ -592,6 +596,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
const struct engine *e = r->e;
+ const float hydro_h_max = e->hydro_properties->h_max;
const float target_wcount = e->hydro_properties->target_neighbours;
const float max_wcount =
target_wcount + e->hydro_properties->delta_neighbours;
@@ -663,15 +668,23 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, correct then */
p->h += h_corr;
- /* Flag for another round of fun */
- pid[redo] = pid[i];
- redo += 1;
+ /* If below the absolute maximum, try again */
+ if (p->h < hydro_h_max) {
+
+ /* Flag for another round of fun */
+ pid[redo] = pid[i];
+ redo += 1;
- /* Re-initialise everything */
- hydro_init_part(p);
+ /* Re-initialise everything */
+ hydro_init_part(p);
+
+ /* Off we go ! */
+ continue;
+ } else {
- /* Off we go ! */
- continue;
+ /* Ok, this particle is a lost cause... */
+ p->h = hydro_h_max;
+ }
}
/* We now have a particle whose smoothing length has converged */
@@ -1337,13 +1350,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Get a handle on the part. */
struct part *restrict p = &parts[k];
-
- if (part_is_active(p, e)) {
-
- /* First, finish the force loop */
- hydro_end_force(p);
- if (p->gpart != NULL) gravity_end_force(p->gpart, const_G);
- }
+ if (part_is_active(p, e)) hydro_end_force(p);
}
/* Loop over the g-particles in this cell. */
@@ -1351,24 +1358,26 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Get a handle on the gpart. */
struct gpart *restrict gp = &gparts[k];
+ if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
- if (gp->type == swift_type_dark_matter) {
- if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
+#ifdef SWIFT_DEBUG_CHECKS
+ if (e->policy & engine_policy_self_gravity) {
+ gp->mass_interacted += gp->mass;
+ if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
+ error(
+ "g-particle did not interact gravitationally with all other "
+ "particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
+ gp->mass_interacted, e->s->total_mass, gp->mass);
}
+#endif
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
- /* Get a handle on the part. */
+ /* Get a handle on the spart. */
struct spart *restrict sp = &sparts[k];
-
- if (spart_is_active(sp, e)) {
-
- /* First, finish the force loop */
- star_end_force(sp);
- gravity_end_force(sp->gpart, const_G);
- }
+ if (spart_is_active(sp, e)) star_end_force(sp);
}
}
@@ -1617,6 +1626,8 @@ void *runner_main(void *data) {
/* Get the cells. */
struct cell *ci = t->ci;
struct cell *cj = t->cj;
+
+/* Mark the thread we run on */
#ifdef SWIFT_DEBUG_TASKS
t->rid = r->cpuid;
#endif
@@ -1688,8 +1699,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_doself2_force(r, ci);
else if (t->subtype == task_subtype_grav)
- // runner_doself_grav(r, ci, 1);
- ;
+ runner_doself_grav(r, ci, 1);
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else
@@ -1706,8 +1716,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dopair2_force(r, ci, cj);
else if (t->subtype == task_subtype_grav)
- // runner_dopair_grav(r, ci, cj, 1);#
- ;
+ runner_dopair_grav(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -1722,8 +1731,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dosub_self2_force(r, ci, 1);
else if (t->subtype == task_subtype_grav)
- // runner_dosub_grav(r, ci, cj, 1);#
- ;
+ runner_dosub_grav(r, ci, cj, 1);
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else
@@ -1740,8 +1748,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dosub_pair2_force(r, ci, cj, t->flags, 1);
else if (t->subtype == task_subtype_grav)
- // runner_dosub_grav(r, ci, cj, 1);
- ;
+ runner_dosub_grav(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
@@ -1799,17 +1806,17 @@ void *runner_main(void *data) {
break;
#endif
case task_type_grav_mm:
- // runner_do_grav_mm(r, t->ci, 1);
- ;
+ // runner_do_grav_mm(r, t->ci, 1);
+ break;
+ case task_type_grav_down:
+ runner_do_grav_down(r, t->ci);
+ break;
+ case task_type_grav_top_level:
+ // runner_do_grav_top_level(r);
+ break;
+ case task_type_grav_long_range:
+ runner_do_grav_long_range(r, t->ci, 1);
break;
- /* case task_type_grav_up: */
- /* runner_do_grav_up(r, t->ci); */
- /* break; */
- /* case task_type_grav_gather_m: */
- /* break; */
- /* case task_type_grav_fft: */
- /* runner_do_grav_fft(r); */
- /* break; */
case task_type_cooling:
if (e->policy & engine_policy_cooling) runner_do_end_force(r, ci, 1);
runner_do_cooling(r, t->ci, 1);
@@ -1821,6 +1828,18 @@ void *runner_main(void *data) {
error("Unknown/invalid task type (%d).", t->type);
}
+/* Mark that we have run this task on these cells */
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci != NULL) {
+ ci->tasks_executed[t->type]++;
+ ci->subtasks_executed[t->subtype]++;
+ }
+ if (cj != NULL) {
+ cj->tasks_executed[t->type]++;
+ cj->subtasks_executed[t->subtype]++;
+ }
+#endif
+
/* We're done with this task, see if we get a next one. */
prev = t;
t = scheduler_done(sched, t);
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9d2606ceb06fd6d32592010376e867a6ae582bf0..9988b7d553a962ee541f20cab64cc534c991957a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -34,25 +34,78 @@
*/
void runner_do_grav_up(struct runner *r, struct cell *c) {
- if (c->split) { /* Regular node */
+ /* if (c->split) { /\* Regular node *\/ */
+
+ /* /\* Recurse. *\/ */
+ /* for (int k = 0; k < 8; k++) */
+ /* if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]); */
+
+ /* /\* Collect the multipoles from the progeny. *\/ */
+ /* multipole_reset(&c->multipole); */
+ /* for (int k = 0; k < 8; k++) { */
+ /* if (c->progeny[k] != NULL) */
+ /* multipole_add(&c->multipole, &c->progeny[k]->multipole); */
+ /* } */
+
+ /* } else { /\* Leaf node. *\/ */
+ /* /\* Just construct the multipole from the gparts. *\/ */
+ /* multipole_init(&c->multipole, c->gparts, c->gcount); */
+ /* } */
+}
+
+void runner_do_grav_down(struct runner *r, struct cell *c) {
+
+ if (c->split) {
- /* Recurse. */
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]);
+ for (int k = 0; k < 8; ++k) {
+ struct cell *cp = c->progeny[k];
+ struct gravity_tensors temp;
- /* Collect the multipoles from the progeny. */
- multipole_reset(&c->multipole);
- for (int k = 0; k < 8; k++) {
- if (c->progeny[k] != NULL)
- multipole_add(&c->multipole, &c->progeny[k]->multipole);
+ if (cp != NULL) {
+ gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
+ 1);
+ }
}
- } else { /* Leaf node. */
- /* Just construct the multipole from the gparts. */
- multipole_init(&c->multipole, c->gparts, c->gcount);
+ } else {
+
+ gravity_L2P(c->multipole, c->gparts, c->gcount);
}
}
+/**
+ * @brief Computes the interaction of the field tensor in a cell with the
+ * multipole of another cell.
+ *
+ * @param r The #runner.
+ * @param ci The #cell with field tensor to interact.
+ * @param cj The #cell with the multipole.
+ */
+__attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
+ const struct runner *r, const struct cell *restrict ci,
+ const struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+ const int periodic = e->s->periodic;
+ const struct multipole *multi_j = &cj->multipole->m_pole;
+ // const float a_smooth = e->gravity_properties->a_smooth;
+ // const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (multi_j->mass == 0.0) error("Multipole does not seem to have been set.");
+#endif
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e)) return;
+
+ gravity_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
+ periodic);
+
+ TIMER_TOC(timer_dopair_grav_mm);
+}
+
/**
* @brief Computes the interaction of all the particles in a cell with the
* multipole of another cell.
@@ -68,15 +121,16 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const struct engine *e = r->e;
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
- const struct multipole multi = cj->multipole;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const struct gravity_tensors *multi = cj->multipole;
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (gcount == 0) error("Empty cell!"); // MATTHIEU sanity check
+ if (gcount == 0) error("Empty cell!");
- if (multi.mass == 0.0) // MATTHIEU sanity check
+ if (multi->m_pole.mass == 0.0)
error("Multipole does not seem to have been set.");
#endif
@@ -105,13 +159,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
if (!gpart_is_active(gp, e)) continue;
/* Compute the pairwise distance. */
- const float dx[3] = {multi.CoM[0] - gp->x[0], // x
- multi.CoM[1] - gp->x[1], // y
- multi.CoM[2] - gp->x[2]}; // z
+ const float dx[3] = {multi->CoM[0] - gp->x[0], // x
+ multi->CoM[1] - gp->x[1], // y
+ multi->CoM[2] - gp->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact !*/
- runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi);
+ runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->mass_interacted += multi->m_pole.mass;
+#endif
}
TIMER_TOC(timer_dopair_grav_pm);
@@ -135,7 +193,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -195,18 +254,30 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+ gpj->mass_interacted += gpi->mass;
+#endif
+
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+#endif
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->mass_interacted += gpi->mass;
+#endif
}
}
}
@@ -229,7 +300,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
TIMER_TIC;
@@ -277,18 +349,31 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+ gpj->mass_interacted += gpi->mass;
+#endif
+
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+#endif
+
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->mass_interacted += gpi->mass;
+#endif
}
}
}
@@ -466,7 +551,8 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
}
}
-static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
+static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+ int timer) {
#if ICHECK > 0
for (int pid = 0; pid < ci->gcount; pid++) {
@@ -485,10 +571,10 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
const struct engine *e = r->e;
struct cell *cells = e->s->cells_top;
const int nr_cells = e->s->nr_cells;
- const double max_d =
- const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
- const double max_d2 = max_d * max_d;
- const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
+ /* const double max_d = */
+ /* const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */
+ /* const double max_d2 = max_d * max_d; */
+ // const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
@@ -501,14 +587,14 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
if (ci == cj) continue;
- const double dx[3] = {cj->loc[0] - pos_i[0], // x
- cj->loc[1] - pos_i[1], // y
- cj->loc[2] - pos_i[2]}; // z
- const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ /* const double dx[3] = {cj->loc[0] - pos_i[0], // x */
+ /* cj->loc[1] - pos_i[1], // y */
+ /* cj->loc[2] - pos_i[2]}; // z */
+ /* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
- if (r2 > max_d2) continue;
+ // if (r2 > max_d2) continue;
- if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
+ if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}
}
diff --git a/src/space.c b/src/space.c
index 6dba35bf60da20d8c545f6c3317ed8de957a643f..625fe944c488c871d02b14c72d15051b3e53b3c4 100644
--- a/src/space.c
+++ b/src/space.c
@@ -186,8 +186,8 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
struct cell **cell_rec_begin,
struct cell **cell_rec_end,
- struct multipole **multipole_rec_begin,
- struct multipole **multipole_rec_end) {
+ struct gravity_tensors **multipole_rec_begin,
+ struct gravity_tensors **multipole_rec_end) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
@@ -221,7 +221,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
for (int k = 0; k < num_elements; k++) {
struct cell *c = &cells[k];
struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
- struct multipole *multipole_rec_begin = NULL, *multipole_rec_end = NULL;
+ struct gravity_tensors *multipole_rec_begin = NULL,
+ *multipole_rec_end = NULL;
space_rebuild_recycle_rec(s, c, &cell_rec_begin, &cell_rec_end,
&multipole_rec_begin, &multipole_rec_end);
if (cell_rec_begin != NULL)
@@ -247,6 +248,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->drift = NULL;
c->cooling = NULL;
c->sourceterms = NULL;
+ c->grav_top_level = NULL;
+ c->grav_long_range = NULL;
+ c->grav_down = NULL;
c->super = c;
if (c->sort != NULL) {
free(c->sort);
@@ -409,9 +413,9 @@ void space_regrid(struct space *s, int verbose) {
/* Allocate the multipoles for the top-level cells. */
if (s->gravity) {
if (posix_memalign((void *)&s->multipoles_top, multipole_align,
- s->nr_cells * sizeof(struct multipole)) != 0)
+ s->nr_cells * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate top-level multipoles.");
- bzero(s->multipoles_top, s->nr_cells * sizeof(struct multipole));
+ bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
}
/* Set the cells' locks */
@@ -2087,7 +2091,7 @@ void space_split_recursive(struct space *s, struct cell *c,
if (s->gravity) {
/* Reset everything */
- multipole_init(c->multipole);
+ gravity_reset(c->multipole);
/* Compute CoM of all progenies */
double CoM[3] = {0., 0., 0.};
@@ -2095,11 +2099,11 @@ void space_split_recursive(struct space *s, struct cell *c,
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
- const struct multipole *m = c->progeny[k]->multipole;
- CoM[0] += m->CoM[0] * m->mass;
- CoM[1] += m->CoM[1] * m->mass;
- CoM[2] += m->CoM[2] * m->mass;
- mass += m->mass;
+ const struct gravity_tensors *m = c->progeny[k]->multipole;
+ CoM[0] += m->CoM[0] * m->m_pole.mass;
+ CoM[1] += m->CoM[1] * m->m_pole.mass;
+ CoM[2] += m->CoM[2] * m->m_pole.mass;
+ mass += m->m_pole.mass;
}
}
c->multipole->CoM[0] = CoM[0] / mass;
@@ -2110,9 +2114,11 @@ void space_split_recursive(struct space *s, struct cell *c,
struct multipole temp;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
- const struct multipole *m = c->progeny[k]->multipole;
- multipole_M2M(&temp, m, c->multipole->CoM, m->CoM, s->periodic);
- multipole_add(c->multipole, &temp);
+ const struct cell *cp = c->progeny[k];
+ const struct multipole *m = &cp->multipole->m_pole;
+ gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
+ s->periodic);
+ gravity_multipole_add(&c->multipole->m_pole, &temp);
}
}
}
@@ -2160,12 +2166,15 @@ void space_split_recursive(struct space *s, struct cell *c,
struct spart *sp = &sparts[k];
const integertime_t ti_end =
get_integer_time_end(e->ti_current, sp->time_bin);
+ const integertime_t ti_beg =
+ get_integer_time_begin(e->ti_current + 1, sp->time_bin);
if (ti_end < ti_end_min) ti_end_min = ti_end;
if (ti_end > ti_end_max) ti_end_max = ti_end;
+ if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
}
/* Construct the multipole and the centre of mass*/
- if (s->gravity) multipole_P2M(c->multipole, c->gparts, c->gcount);
+ if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
}
/* Set the values for this cell. */
@@ -2278,8 +2287,8 @@ void space_recycle(struct space *s, struct cell *c) {
*/
void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct cell *cell_list_end,
- struct multipole *multipole_list_begin,
- struct multipole *multipole_list_end) {
+ struct gravity_tensors *multipole_list_begin,
+ struct gravity_tensors *multipole_list_end) {
int count = 0;
@@ -2348,8 +2357,9 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Is the multipole buffer empty? */
if (s->gravity && s->multipoles_sub == NULL) {
- if (posix_memalign((void *)&s->multipoles_sub, multipole_align,
- space_cellallocchunk * sizeof(struct multipole)) != 0)
+ if (posix_memalign(
+ (void *)&s->multipoles_sub, multipole_align,
+ space_cellallocchunk * sizeof(struct gravity_tensors)) != 0)
error("Failed to allocate more multipoles.");
/* Constructed a linked list */
@@ -2375,7 +2385,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Init some things in the cell we just got. */
for (int j = 0; j < nr_cells; j++) {
- struct multipole *temp = cells[j]->multipole;
+ struct gravity_tensors *temp = cells[j]->multipole;
bzero(cells[j], sizeof(struct cell));
cells[j]->multipole = temp;
cells[j]->nodeID = -1;
@@ -2827,27 +2837,51 @@ void space_link_cleanup(struct space *s) {
/**
* @brief Checks that all cells have been drifted to a given point in time
*
- * Expensive function. Should only be used for debugging purposes.
+ * Should only be used for debugging purposes.
*
* @param s The #space to check.
* @param ti_drift The (integer) time.
*/
void space_check_drift_point(struct space *s, integertime_t ti_drift) {
-
+#ifdef SWIFT_DEBUG_CHECKS
/* Recursively check all cells */
space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
}
/**
* @brief Checks that all particles and local cells have a non-zero time-step.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to check.
*/
void space_check_timesteps(struct space *s) {
#ifdef SWIFT_DEBUG_CHECKS
-
for (int i = 0; i < s->nr_cells; ++i) {
cell_check_timesteps(&s->cells_top[i]);
}
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+/**
+ * @brief Resets all the individual cell task counters to 0.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to reset.
+ */
+void space_reset_task_counters(struct space *s) {
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int i = 0; i < s->nr_cells; ++i) {
+ cell_reset_task_counters(&s->cells_top[i]);
+ }
+#else
+ error("Calling debugging code without debugging flag activated.");
#endif
}
diff --git a/src/space.h b/src/space.h
index 4bbb8f32353c795772a461061096a366d0848f7e..e3886364bfa7b2f85e724993e5fe3c2d30663bd2 100644
--- a/src/space.h
+++ b/src/space.h
@@ -72,6 +72,9 @@ struct space {
/*! Are we doing gravity? */
int gravity;
+ /*! Total mass in the system */
+ double total_mass;
+
/*! Width of the top-level cells. */
double width[3];
@@ -103,10 +106,10 @@ struct space {
struct cell *cells_sub;
/*! The multipoles associated with the top-level (level 0) cells */
- struct multipole *multipoles_top;
+ struct gravity_tensors *multipoles_top;
/*! Buffer of unused multipoles for the sub-cells. */
- struct multipole *multipoles_sub;
+ struct gravity_tensors *multipoles_sub;
/*! The total number of parts in the space. */
size_t nr_parts, size_parts;
@@ -194,8 +197,8 @@ void space_rebuild(struct space *s, int verbose);
void space_recycle(struct space *s, struct cell *c);
void space_recycle_list(struct space *s, struct cell *cell_list_begin,
struct cell *cell_list_end,
- struct multipole *multipole_list_begin,
- struct multipole *multipole_list_end);
+ struct gravity_tensors *multipole_list_begin,
+ struct gravity_tensors *multipole_list_end);
void space_split(struct space *s, struct cell *cells, int nr_cells,
int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
@@ -215,6 +218,7 @@ void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
+void space_reset_task_counters(struct space *s);
void space_clean(struct space *s);
#endif /* SWIFT_SPACE_H */
diff --git a/src/timeline.h b/src/timeline.h
index c73b2432b219a8ab0254d21c59102841557a57b9..25ab231e24b85f2a4b6b3d44f024588420e432ea 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -37,6 +37,8 @@ typedef char timebin_t;
/*! The maximal number of timesteps in a simulation */
#define max_nr_timesteps (1LL << (num_time_bins + 1))
+#define time_bin_inactive (num_time_bins + 2)
+
/**
* @brief Returns the integer time interval corresponding to a time bin
*
diff --git a/src/timers.h b/src/timers.h
index 50f630e7fc355808596967d8d7d887583674d24a..4cb4d7e0ba60003ba4caefffe257c929c59a8d9e 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -46,6 +46,7 @@ enum {
timer_dopair_gradient,
timer_dopair_force,
timer_dopair_grav_pm,
+ timer_dopair_grav_mm,
timer_dopair_grav_pp,
timer_dograv_external,
timer_dosource,
diff --git a/src/timestep.h b/src/timestep.h
index 6497e25984ef36d759afd12616b45c8166816dfb..99ca1b4f90cc7b8894d4dfe0e0d2670da8f01d65 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -78,7 +78,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
e->physical_constants, gp));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(gp));
+ new_dt =
+ min(new_dt, gravity_compute_timestep_self(gp, e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
@@ -121,7 +122,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
e->physical_constants, p->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
+ new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(
+ p->gpart, e->gravity_properties));
}
/* Final time-step is minimum of hydro and gravity */
@@ -163,7 +165,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
e->physical_constants, sp->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart));
+ new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart,
+ e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
diff --git a/src/tools.c b/src/tools.c
index ab11d1f5930cf5319aaf6424f1559f144718e154..89ac286fb435c01b361bdea66e62dd2d7f41ee24 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -732,19 +732,24 @@ int compare_particles(struct part a, struct part b, double threshold) {
#endif
}
-/** @brief Computes the forces between all g-particles using the N^2 algorithm
+/**
+ * @brief Computes the forces between all g-particles using the N^2 algorithm
*
* Overwrites the accelerations of the gparts with the values.
* Do not use for actual runs.
*
* @brief gparts The array of particles.
* @brief gcount The number of particles.
+ * @brief constants Physical constants in internal units.
+ * @brief gravity_properties Constants governing the gravity scheme.
*/
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr) {
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
- const float max_d = const_gravity_r_cut * rlr;
+ const float r_cut = gravity_properties->r_cut;
+ const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
message("rlr_inv= %f", rlr_inv);
diff --git a/src/tools.h b/src/tools.h
index ece3078dce7cc8ab4b15538a1e5d9a990d81b36d..4d9e8d3ef86f9ad2661118acf008797893ea5bd7 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -23,6 +23,7 @@
#define SWIFT_TOOL_H
#include "cell.h"
+#include "gravity_properties.h"
#include "part.h"
#include "physical_constants.h"
#include "runner.h"
@@ -45,8 +46,8 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
double random_uniform(double a, double b);
void shuffle_particles(struct part *parts, const int count);
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr);
-
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr);
int compare_values(double a, double b, double threshold, double *absDiff,
double *absSum, double *relDiff);
int compare_particles(struct part a, struct part b, double threshold);
diff --git a/tests/testKernelGrav.c b/tests/testKernelGrav.c
index 2733a4a4d041149499d354ef217ae85b1cd35b7f..41e085945693efaf658c219d0e5f992ddf023d74 100644
--- a/tests/testKernelGrav.c
+++ b/tests/testKernelGrav.c
@@ -87,7 +87,7 @@ int main() {
printf("\nAll values are consistent\n");
/* Now test the long range function */
- const float a_smooth = const_gravity_a_smooth;
+ const float a_smooth = 4.5f;
for (int k = 1; k < numPoints; ++k) {