Commit d44fbb02 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' into initialisation_fixes

parents 78b42529 ad5c2154
......@@ -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
*
......@@ -2305,7 +2440,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 */
......@@ -2316,6 +2460,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;
......@@ -2440,7 +2585,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 */
......@@ -2450,6 +2604,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
*
* This function only branches towards the potential chosen by the user.