Commit ad5c2154 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'gravity_infrastructure' into 'master'

Gravity infrastructure

Sorry, this is a bulky one. Not all of the gravity is here but most of the infrastructure is in place. It's getting tricky to maintain this separate so I'd like to merge all the innocuous changes to the master branch.
All of the gravity tasks and the logic is in place. Only the content of some of them needs to be more thoroughly checked for accuracy and I need to merge the FFT task once I am happy with it. 

With this in it should be much easier for me and Bert to re-import gizmo into the master.

Could you stress test this and verify that your usual hydro test-cases still run ? I'll bring gravity test cases in at a later point.

See merge request !212
parents 6f26c406 3411dd0f
......@@ -56,6 +56,8 @@ tests/testSingle
tests/testTimeIntegration
tests/testSPHStep
tests/testKernel
tests/testKernelGrav
tests/testFFT
tests/testInteractions
tests/testSymmetry
tests/testMaths
......
......@@ -417,6 +417,17 @@ if test "$ac_cv_func_pthread_setaffinity_np" = "yes"; then
fi
fi
# Check for FFTW
have_fftw3="no"
AC_CHECK_HEADER([fftw3.h])
if test "$ac_cv_header_fftw3_h" = "yes"; then
AC_CHECK_LIB([fftw3],[fftw_malloc], [AC_DEFINE([HAVE_FFTW],
[1],[Defined if FFTW 3.x exists.])] )
FFTW_LIBS="-lfftw3"
have_fftw3="yes"
fi
AC_SUBST([FFTW_LIBS])
# Check for Intel intrinsics header optionally used by vector.h.
AC_CHECK_HEADERS([immintrin.h])
......@@ -503,6 +514,7 @@ AC_MSG_RESULT([
HDF5 enabled : $with_hdf5
- parallel : $have_parallel_hdf5
Metis enabled : $have_metis
FFTW3 enabled : $have_fftw3
libNUMA enabled : $have_numa
Using tcmalloc : $have_tcmalloc
CPU profiler : $have_profiler
......
......@@ -24,7 +24,7 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
AM_LDFLAGS = $(HDF5_LDFLAGS)
# Extra libraries.
EXTRA_LIBS = $(HDF5_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS)
EXTRA_LIBS = $(HDF5_LIBS) $(FFTW_LIBS) $(PROFILER_LIBS) $(TCMALLOC_LIBS)
# MPI libraries.
MPI_LIBS = $(METIS_LIBS) $(MPI_THREAD_LIBS)
......
......@@ -26,7 +26,7 @@ from numpy import *
# with a density of 1
# Parameters
periodic= 1 # 1 For periodic box
periodic= 0 # 1 For periodic box
boxSize = 1.
rho = 1.
L = int(sys.argv[1]) # Number of particles along one axis
......
......@@ -399,8 +399,8 @@ int main(int argc, char *argv[]) {
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
struct space s;
space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic, talking,
dry_run);
space_init(&s, params, dim, parts, gparts, Ngas, Ngpart, periodic,
with_self_gravity, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
......
......@@ -19,7 +19,7 @@
AM_CFLAGS = -DTIMER $(HDF5_CPPFLAGS)
# Assign a "safe" version number
AM_LDFLAGS = $(HDF5_LDFLAGS) -version-info 0:0:0
AM_LDFLAGS = $(HDF5_LDFLAGS) $(FFTW_LIBS) -version-info 0:0:0
# The git command, if available.
GIT_CMD = @GIT_CMD@
......@@ -48,17 +48,18 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potentials.c hydro_properties.c
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potentials.c hydro_properties.c \
runner_doiact_fft.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \
timestep.h drift.h adiabatic_index.h io_properties.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
hydro.h hydro_io.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
hydro/Default/hydro.h hydro/Default/hydro_iact.h hydro/Default/hydro_io.h \
......
......@@ -709,6 +709,92 @@ void cell_clean_links(struct cell *c, void *data) {
}
/**
<<<<<<< HEAD
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
*
* @param ci First #cell.
* @param cj Second #cell.
*
* @todo Deal with periodicity.
*/
int cell_are_neighbours(const struct cell *restrict ci,
const struct cell *restrict cj) {
#ifdef SWIFT_DEBUG_CHECKS
if (ci->width[0] != cj->width[0]) error("Cells of different size !");
#endif
/* Maximum allowed distance */
const double min_dist =
1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
/* (Manhattan) Distance between the cells */
for (int k = 0; k < 3; k++) {
const double center_i = ci->loc[k];
const double center_j = cj->loc[k];
if (fabsf(center_i - center_j) > min_dist) return 0;
}
return 1;
}
/**
* @brief Computes the multi-pole brutally and compare to the
* recursively computed one.
*
* @param c Cell to act upon
* @param data Unused parameter
*/
void cell_check_multipole(struct cell *c, void *data) {
struct multipole ma;
if (c->gcount > 0) {
/* Brute-force calculation */
multipole_init(&ma, c->gparts, c->gcount);
/* Compare with recursive one */
struct multipole mb = c->multipole;
if (fabsf(ma.mass - mb.mass) / fabsf(ma.mass + mb.mass) > 1e-5)
error("Multipole masses are different (%12.15e vs. %12.15e)", ma.mass,
mb.mass);
for (int k = 0; k < 3; ++k)
if (fabsf(ma.CoM[k] - mb.CoM[k]) / fabsf(ma.CoM[k] + mb.CoM[k]) > 1e-5)
error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
mb.CoM[k]);
if (fabsf(ma.I_xx - mb.I_xx) / fabsf(ma.I_xx + mb.I_xx) > 1e-5 &&
ma.I_xx > 1e-9)
error("Multipole I_xx are different (%12.15e vs. %12.15e)", ma.I_xx,
mb.I_xx);
if (fabsf(ma.I_yy - mb.I_yy) / fabsf(ma.I_yy + mb.I_yy) > 1e-5 &&
ma.I_yy > 1e-9)
error("Multipole I_yy are different (%12.15e vs. %12.15e)", ma.I_yy,
mb.I_yy);
if (fabsf(ma.I_zz - mb.I_zz) / fabsf(ma.I_zz + mb.I_zz) > 1e-5 &&
ma.I_zz > 1e-9)
error("Multipole I_zz are different (%12.15e vs. %12.15e)", ma.I_zz,
mb.I_zz);
if (fabsf(ma.I_xy - mb.I_xy) / fabsf(ma.I_xy + mb.I_xy) > 1e-5 &&
ma.I_xy > 1e-9)
error("Multipole I_xy are different (%12.15e vs. %12.15e)", ma.I_xy,
mb.I_xy);
if (fabsf(ma.I_xz - mb.I_xz) / fabsf(ma.I_xz + mb.I_xz) > 1e-5 &&
ma.I_xz > 1e-9)
error("Multipole I_xz are different (%12.15e vs. %12.15e)", ma.I_xz,
mb.I_xz);
if (fabsf(ma.I_yz - mb.I_yz) / fabsf(ma.I_yz + mb.I_yz) > 1e-5 &&
ma.I_yz > 1e-9)
error("Multipole I_yz are different (%12.15e vs. %12.15e)", ma.I_yz,
mb.I_yz);
}
}
/*
* @brief Frees up the memory allocated for this #cell
*/
void cell_clean(struct cell *c) {
......
......@@ -197,6 +197,9 @@ void cell_init_parts(struct cell *c, void *data);
void cell_init_gparts(struct cell *c, void *data);
void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell *c, void *data);
int cell_are_neighbours(const struct cell *restrict ci,
const struct cell *restrict cj);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
#endif /* SWIFT_CELL_H */
......@@ -374,6 +374,9 @@ void writeCodeDescription(hid_t h_file) {
writeAttribute_s(h_grpcode, "Git Branch", git_branch());
writeAttribute_s(h_grpcode, "Git Revision", git_revision());
writeAttribute_s(h_grpcode, "HDF5 library version", hdf5_version());
#ifdef HAVE_FFTW
writeAttribute_s(h_grpcode, "FFTW library version", fftw3_version());
#endif
#ifdef WITH_MPI
writeAttribute_s(h_grpcode, "MPI library", mpi_version());
#ifdef HAVE_METIS
......
......@@ -36,11 +36,6 @@
/* Time integration constants. */
#define const_max_u_change 0.1f
/* Gravity stuff. */
#define const_theta_max 0.57735f
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
/* Hydrodynamical adiabatic index. */
#define HYDRO_GAMMA_5_3
//#define HYDRO_GAMMA_4_3
......@@ -59,7 +54,13 @@
#define GADGET2_SPH
//#define DEFAULT_SPH
/* External potential properties */
/* Self gravity stuff. */
#define const_gravity_multipole_order 2
#define const_gravity_a_smooth 1.25f
#define const_gravity_r_cut 4.5f
#define const_gravity_eta 0.025f
/* External gravity properties */
#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
......
......@@ -26,6 +26,7 @@
#include "const.h"
#include "debug.h"
#include "hydro.h"
#include "part.h"
/**
* @brief Perform the 'drift' operation on a #gpart
......
......@@ -64,6 +64,7 @@
#include "serial_io.h"
#include "single_io.h"
#include "timers.h"
#include "tools.h"
#include "units.h"
#include "version.h"
......@@ -129,12 +130,10 @@ void engine_make_gravity_hierarchical_tasks(struct engine *e, struct cell *c,
engine_policy_external_gravity;
const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
/* Am I the super-cell? */
/* TODO(pedro): Add a condition for gravity tasks as well. */
if (super == NULL &&
(c->density != NULL || (!c->split && (c->count > 0 || c->gcount > 0)))) {
/* Is this the super-cell? */
if (super == NULL && (c->grav != NULL || (!c->split && c->gcount > 0))) {
/* This is the super cell, i.e. the first with density tasks attached. */
/* This is the super cell, i.e. the first with gravity tasks attached. */
super = c;
/* Local tasks only... */
......@@ -1182,9 +1181,62 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
#endif
}
/**
* @brief Constructs the top-level tasks for the short-range gravity
* interactions.
*
* All top-cells get a self task.
* All neighbouring pairs get a pair task.
* All non-neighbouring pairs within a range of 6 cells get a M-M task.
*
* @param e The #engine.
*/
void engine_make_gravity_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
struct cell *cells = s->cells;
const int nr_cells = s->nr_cells;
for (int cid = 0; cid < nr_cells; ++cid) {
struct cell *ci = &cells[cid];
/* Skip cells without gravity particles */
if (ci->gcount == 0) continue;
/* Is that neighbour local ? */
if (ci->nodeID != nodeID) continue;
/* If the cells is local build a self-interaction */
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];
/* Skip cells without gravity particles */
if (cj->gcount == 0) continue;
/* Is that neighbour local ? */
if (cj->nodeID != nodeID) continue;
if (cell_are_neighbours(ci, cj))
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
cj, 1);
}
}
}
/**
* @brief Constructs the top-level pair tasks for the first hydro loop over
*neighbours
* neighbours
*
* Here we construct all the tasks for all possible neighbouring non-empty
* local cells in the hierarchy. No dependencies are being added thus far.
......@@ -1316,6 +1368,115 @@ 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->kick);
/* grav_up --> gravity ( --> kick) */
scheduler_addunlock(sched, c->super->grav_up, gravity);
}
/**
* @brief Creates all the task dependencies for the gravity
*
* @param e The #engine
*/
void engine_link_gravity_tasks(struct engine *e) {
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int nr_tasks = sched->nr_tasks;
/* Add one task gathering all the multipoles */
struct task *gather = scheduler_addtask(
sched, task_type_grav_gather_m, task_subtype_none, 0, 0, NULL, NULL, 0);
/* And one task performing the FFT */
struct task *fft = scheduler_addtask(sched, task_type_grav_fft,
task_subtype_none, 0, 0, NULL, NULL, 0);
scheduler_addunlock(sched, gather, fft);
for (int k = 0; k < nr_tasks; k++) {
/* Get a pointer to the task. */
struct task *t = &sched->tasks[k];
/* Skip? */
if (t->skip) continue;
/* Multipole construction */
if (t->type == task_type_grav_up) {
scheduler_addunlock(sched, t, gather);
}
/* Long-range interaction */
if (t->type == task_type_grav_mm) {
/* Gather the multipoles --> mm interaction --> kick */
scheduler_addunlock(sched, gather, t);
scheduler_addunlock(sched, t, t->ci->super->kick);
/* init --> mm interaction */
scheduler_addunlock(sched, t->ci->super->init, t);
}
/* Self-interaction? */
if (t->type == task_type_self && t->subtype == task_subtype_grav) {
engine_make_gravity_dependencies(sched, t, t->ci);
}
/* Otherwise, pair interaction? */
else if (t->type == task_type_pair && t->subtype == task_subtype_grav) {
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_gravity_dependencies(sched, t, t->cj);
}
}
/* 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);
}
}
/* Otherwise, sub-pair interaction? */
else if (t->type == task_type_sub_pair && t->subtype == task_subtype_grav) {
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_gravity_dependencies(sched, t, t->cj);
}
}
}
}
/**
* @brief Creates the dependency network for the hydro tasks of a given cell.
*
......@@ -1453,41 +1614,6 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
}
}
/**
* @brief Constructs the top-level pair tasks for the gravity M-M interactions
*
* Correct implementation is still lacking here.
*
* @param e The #engine.
*/
void engine_make_gravityinteraction_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
const int nr_cells = s->nr_cells;
struct cell *cells = s->cells;
/* Loop over all cells. */
for (int i = 0; i < nr_cells; i++) {
/* If it has gravity particles, add a self-task */
if (cells[i].gcount > 0) {
scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
&cells[i], NULL, 0);
/* Loop over all remainding cells */
for (int j = i + 1; j < nr_cells; j++) {
/* If that other cell has gravity parts, add a pair interaction */
if (cells[j].gcount > 0) {
scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
&cells[i], &cells[j], 0);
}
}
}
}
}
/**
* @brief Constructs the gravity tasks building the multipoles and propagating
*them to the children
......@@ -1513,9 +1639,12 @@ void engine_make_gravityrecursive_tasks(struct engine *e) {
struct task *up =
scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0, 0,
&cells[k], NULL, 0);
struct task *down =
scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0, 0,
&cells[k], NULL, 0);
struct task *down = NULL;
/* struct task *down = */
/* 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);
......@@ -1548,11 +1677,10 @@ void engine_maketasks(struct engine *e) {
}
/* Construct the firt hydro loop over neighbours */
engine_make_hydroloop_tasks(e);
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_gravityinteraction_tasks(e);
if (e->policy & engine_policy_self_gravity) engine_make_gravity_tasks(e);
/* Split the tasks. */
scheduler_splittasks(sched);
......@@ -1573,7 +1701,7 @@ void engine_maketasks(struct engine *e) {
/* Count the number of tasks associated with each cell and
store the density tasks in each cell, and make each sort
depend on the sorts of its progeny. */
engine_count_and_link_tasks(e);
if (e->policy & engine_policy_hydro) engine_count_and_link_tasks(e);
/* Append hierarchical tasks to each cells */
if (e->policy & engine_policy_hydro)
......@@ -1588,7 +1716,12 @@ void engine_maketasks(struct engine *e) {
/* 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
of its super-cell. */
engine_make_extra_hydroloop_tasks(e);
if (e->policy & engine_policy_hydro) engine_make_extra_hydroloop_tasks(e);
/* Add the dependencies for the self-gravity stuff */
if (e->policy & engine_policy_self_gravity) engine_link_gravity_tasks(e);
#ifdef WITH_MPI
/* Add the communication tasks if MPI is being used. */
if (e->policy & engine_policy_mpi) {
......@@ -1611,6 +1744,7 @@ void engine_maketasks(struct engine *e) {
NULL);
}
}
#endif
/* Set the unlocks per task. */
scheduler_set_unlocks(sched);
......@@ -1740,7 +1874,7 @@ int engine_marktasks(struct engine *e) {
continue;
/* Set the sort flags. */
if (t->type == task_type_pair) {
if (t->type == task_type_pair && t->subtype != task_subtype_grav) {
if (!(ci->sorted & (1 << t->flags))) {
ci->sorts->flags |= (1 << t->flags);
ci->sorts->skip = 0;
......@@ -2144,6 +2278,7 @@ void engine_collect_drift(struct cell *c) {
c->ang_mom[1] = ang_mom[1];
c->ang_mom[2] = ang_mom[2];
}
/**
* @brief Print the conserved quantities statistics to a log file
*
......@@ -2315,7 +2450,16 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Add the tasks corresponding to self-gravity to the masks */
if (e->policy & engine_policy_self_gravity) {
/* Nothing here for now */
mask |= 1 << task_type_grav_up;
mask |= 1 << task_type_grav_mm;
mask |= 1 << task_type_grav_gather_m;
mask |= 1 << task_type_grav_fft;
mask |= 1 << task_type_self;
mask |= 1 << task_type_pair;
mask |= 1 << task_type_sub_self;
mask |= 1 << task_type_sub_pair;
submask |= 1 << task_subtype_grav;
}
/* Add the tasks corresponding to external gravity to the masks */
......@@ -2326,6 +2470,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
/* Add MPI tasks if need be */
if (e->policy & engine_policy_mpi) {
mask |= 1 << task_type_send;
mask |= 1 << task_type_recv;
submask |= 1 << task_subtype_tend;
......@@ -2450,7 +2595,16 @@ void engine_step(struct engine *e) {
/* Add the tasks corresponding to self-gravity to the masks */
if (e->policy & engine_policy_self_gravity) {
/* Nothing here for now */
mask |= 1 << task_type_grav_up;
mask |= 1 << task_type_grav_mm;
mask |= 1 << task_type_grav_gather_m;
mask |= 1 << task_type_grav_fft;
mask |= 1 << task_type_self;
mask |= 1 << task_type_pair;
mask |= 1 << task_type_sub_self;
mask |= 1 << task_type_sub_pair;
submask |= 1 << task_subtype_grav;
}
/* Add the tasks corresponding to external gravity to the masks */
......@@ -2460,6 +2614,7 @@ void engine_step(struct engine *e) {
/* Add MPI tasks if need be */
if (e->policy & engine_policy_mpi) {
mask |= 1 << task_type_send;
mask |= 1 << task_type_recv;
submask |= 1 << task_subtype_tend;
......
......@@ -53,8 +53,6 @@ gravity_compute_timestep_external(const struct external_potential* potential,
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*