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

Merge branch 'gravity_multi_dt' into 'master'

Gravity multi dt

This implements a small change to the top-level multipole drifting strategy. 

It should have no side effect on the hydro business but I'd rather be sure. Could you check that all your usual cases are not broken by this ? Thanks !

See merge request !328
parents 80610af7 08b2c5fb
......@@ -88,6 +88,8 @@ theory/SPH/Flavours/sph_flavours.pdf
theory/SPH/EoS/eos.pdf
theory/SPH/*.pdf
theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf
m4/libtool.m4
m4/ltoptions.m4
......
......@@ -1017,14 +1017,15 @@ void cell_clean_links(struct cell *c, void *data) {
}
/**
* @brief Checks that a cell is at the current point in time
* @brief Checks that the particles in a cell are at the
* current point in time
*
* Calls error() if the cell is not at the current time.
*
* @param c Cell to act upon
* @param data The current time on the integer time-line
*/
void cell_check_drift_point(struct cell *c, void *data) {
void cell_check_particle_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
......@@ -1056,6 +1057,34 @@ void cell_check_drift_point(struct cell *c, void *data) {
#endif
}
/**
* @brief Checks that the multipole of a cell is at the current point in time
*
* Calls error() if the cell is not at the current time.
*
* @param c Cell to act upon
* @param data The current time on the integer time-line
*/
void cell_check_multipole_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_drift = *(integertime_t *)data;
/* Only check local cells */
if (c->nodeID != engine_rank) return;
if (c->ti_old_multipole != ti_drift)
error(
"Cell multipole in an incorrect time-zone! c->ti_old_multipole=%lld "
"ti_drift=%lld (depth=%d)",
c->ti_old_multipole, ti_drift, c->depth);
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
/**
* @brief Resets all the individual cell task counters to 0.
*
......@@ -1128,7 +1157,8 @@ void cell_check_multipole(struct cell *c, void *data) {
/* Now compare the multipole expansion */
if (!gravity_multipole_equal(&ma, c->multipole, tolerance)) {
message("Multipoles are not equal at depth=%d!", c->depth);
message("Multipoles are not equal at depth=%d! tol=%f", c->depth,
tolerance);
message("Correct answer:");
gravity_multipole_print(&ma.m_pole);
message("Recursive multipole:");
......@@ -1476,17 +1506,15 @@ void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
/* Check that we are actually going to move forward. */
if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
/* Drift the multipole */
if (ti_current > ti_old_multipole) gravity_drift(c->multipole, dt);
/* Are we not in a leaf ? */
if (c->split) {
/* Loop over the progeny and drift the multipoles. */
/* Loop over the progeny and recurse. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) cell_drift_particles(c->progeny[k], e);
} else if (ti_current > ti_old_multipole) {
/* Drift the multipole */
gravity_drift(c->multipole, dt);
if (c->progeny[k] != NULL) cell_drift_all_multipoles(c->progeny[k], e);
}
/* Update the time of the last drift */
......
......@@ -352,7 +352,8 @@ 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);
void cell_check_drift_point(struct cell *c, void *data);
void cell_check_particle_drift_point(struct cell *c, void *data);
void cell_check_multipole_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);
......
......@@ -846,7 +846,8 @@ void engine_repartition(struct engine *e) {
fflush(stdout);
/* Check that all cells have been drifted to the current time */
space_check_drift_point(e->s, e->ti_old);
space_check_drift_point(e->s, e->ti_old,
e->policy & engine_policy_self_gravity);
#endif
/* Clear the repartition flag. */
......@@ -2664,7 +2665,8 @@ void engine_rebuild(struct engine *e) {
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_old);
space_check_drift_point(e->s, e->ti_old,
e->policy & engine_policy_self_gravity);
#endif
if (e->verbose)
......@@ -2686,7 +2688,8 @@ void engine_prepare(struct engine *e) {
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_old);
space_check_drift_point(e->s, e->ti_old,
e->policy & engine_policy_self_gravity);
}
#endif
......@@ -2920,7 +2923,8 @@ void engine_print_stats(struct engine *e) {
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_current);
space_check_drift_point(e->s, e->ti_current,
e->policy & engine_policy_self_gravity);
/* Be verbose about this */
message("Saving statistics at t=%e.", e->time);
......@@ -3219,7 +3223,10 @@ void engine_step(struct engine *e) {
for (int i = 0; i < e->s->nr_cells; ++i)
num_gpart_mpole += e->s->cells_top[i].multipole->m_pole.num_gpart;
if (num_gpart_mpole != e->s->nr_gparts)
error("Multipoles don't contain the total number of gpart");
error(
"Multipoles don't contain the total number of gpart mpoles=%zd "
"ngparts=%zd",
num_gpart_mpole, e->s->nr_gparts);
}
#endif
......@@ -3228,6 +3235,9 @@ void engine_step(struct engine *e) {
gravity_exact_force_compute(e->s, e);
#endif
/* Do we need to drift the top-level multipoles ? */
if (e->policy & engine_policy_self_gravity) engine_drift_top_multipoles(e);
/* Start all the tasks. */
TIMER_TIC;
engine_launch(e, e->nr_threads);
......@@ -3238,6 +3248,10 @@ void engine_step(struct engine *e) {
gravity_exact_force_check(e->s, e, 1e-1);
#endif
if (!(e->policy & engine_policy_hydro) &&
(e->policy & engine_policy_self_gravity) && e->step % 20 == 0)
e->forcerebuild = 1;
/* Collect the values of rebuild from all nodes and recover the (integer)
* end of the next time-step. Do these together to reduce the collective MPI
* calls per step, but some of the gathered information is not applied just
......@@ -3317,8 +3331,8 @@ void engine_unskip(struct engine *e) {
}
/**
* @brief Mapper function to drift ALL particle types and multipoles forward in
* time.
* @brief Mapper function to drift *all* particle types and multipoles forward
* in time.
*
* @param map_data An array of #cell%s.
* @param num_elements Chunk size.
......@@ -3336,7 +3350,7 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
/* Drift all the particles */
cell_drift_particles(c, e);
/* Drift the multipole */
/* Drift the multipoles */
if (e->policy & engine_policy_self_gravity)
cell_drift_all_multipoles(c, e);
}
......@@ -3344,7 +3358,8 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
}
/**
* @brief Drift *all* particles forward to the current time.
* @brief Drift *all* particles and multipoles at all levels
* forward to the current time.
*
* @param e The #engine.
*/
......@@ -3361,7 +3376,54 @@ void engine_drift_all(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time. */
space_check_drift_point(e->s, e->ti_current);
space_check_drift_point(e->s, e->ti_current,
e->policy & engine_policy_self_gravity);
#endif
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
/**
* @brief Mapper function to drift *all* top-level multipoles forward in
* time.
*
* @param map_data An array of #cell%s.
* @param num_elements Chunk size.
* @param extra_data Pointer to an #engine.
*/
void engine_do_drift_top_multipoles_mapper(void *map_data, int num_elements,
void *extra_data) {
struct engine *e = (struct engine *)extra_data;
struct cell *cells = (struct cell *)map_data;
for (int ind = 0; ind < num_elements; ind++) {
struct cell *c = &cells[ind];
if (c != NULL && c->nodeID == e->nodeID) {
/* Drift the multipole at this level only */
cell_drift_multipole(c, e);
}
}
}
/**
* @brief Drift *all* top-level multipoles forward to the current time.
*
* @param e The #engine.
*/
void engine_drift_top_multipoles(struct engine *e) {
const ticks tic = getticks();
threadpool_map(&e->threadpool, engine_do_drift_top_multipoles_mapper,
e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that all cells have been drifted to the current time. */
space_check_top_multipoles_drift_point(e->s, e->ti_current);
#endif
if (e->verbose)
......@@ -3577,7 +3639,8 @@ void engine_dump_snapshot(struct engine *e) {
/* Check that all cells have been drifted to the current time.
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_current);
space_check_drift_point(e->s, e->ti_current,
e->policy & engine_policy_self_gravity);
/* Be verbose about this */
message("writing snapshot at t=%e.", e->time);
......
......@@ -255,6 +255,7 @@ void engine_barrier(struct engine *e, int tid);
void engine_compute_next_snapshot_time(struct engine *e);
void engine_unskip(struct engine *e);
void engine_drift_all(struct engine *e);
void engine_drift_top_multipoles(struct engine *e);
void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
......
......@@ -354,6 +354,11 @@ INLINE static void gravity_field_tensors_print(const struct grav_tensor *l) {
#endif
}
INLINE static void gravity_multipole_init(struct multipole *m) {
bzero(m, sizeof(struct multipole));
}
/**
* @brief Prints the content of a #multipole to stdout.
*
......
......@@ -2150,7 +2150,16 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_beg_max = get_integer_time_begin(e->ti_current + 1, time_bin_max);
/* Construct the multipole and the centre of mass*/
if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
if (s->gravity) {
if (gcount > 0)
gravity_P2M(c->multipole, c->gparts, c->gcount);
else {
gravity_multipole_init(&c->multipole->m_pole);
c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
}
}
}
/* Set the values for this cell. */
......@@ -2819,11 +2828,26 @@ void space_link_cleanup(struct space *s) {
*
* @param s The #space to check.
* @param ti_drift The (integer) time.
* @param multipole Are we also checking the multipoles ?
*/
void space_check_drift_point(struct space *s, integertime_t ti_drift) {
void space_check_drift_point(struct space *s, integertime_t ti_drift,
int multipole) {
#ifdef SWIFT_DEBUG_CHECKS
/* Recursively check all cells */
space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
space_map_cells_pre(s, 1, cell_check_particle_drift_point, &ti_drift);
if (multipole)
space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
#else
error("Calling debugging code without debugging flag activated.");
#endif
}
void space_check_top_multipoles_drift_point(struct space *s,
integertime_t ti_drift) {
#ifdef SWIFT_DEBUG_CHECKS
for (int i = 0; i < s->nr_cells; ++i) {
cell_check_multipole_drift_point(&s->cells_top[i], &ti_drift);
}
#else
error("Calling debugging code without debugging flag activated.");
#endif
......
......@@ -213,7 +213,10 @@ void space_init_parts(struct space *s);
void space_init_gparts(struct space *s);
void space_init_sparts(struct space *s);
void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, integertime_t ti_drift);
void space_check_drift_point(struct space *s, integertime_t ti_drift,
int multipole);
void space_check_top_multipoles_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);
......
@ARTICLE{Dehnen2014,
author = {{Dehnen}, W.},
title = "{A fast multipole method for stellar dynamics}",
journal = {Computational Astrophysics and Cosmology},
archivePrefix = "arXiv",
eprint = {1405.2255},
primaryClass = "astro-ph.IM",
year = 2014,
month = sep,
volume = 1,
eid = {1},
pages = {1},
doi = {10.1186/s40668-014-0001-7},
adsurl = {http://adsabs.harvard.edu/abs/2014ComAC...1....1D},
adsnote = {Provided by the SAO/NASA Astrophysics Data System}
}
@article{Hardy2006,
title={Combinatorics of partial derivatives},
author={Hardy, Michael},
journal={Electron. J. Combin},
volume={13},
number={1},
pages={13},
year={2006}
}
\documentclass[fleqn, usenatbib, useAMS,a4paper]{mnras}
\usepackage{graphicx}
\usepackage{amsmath,paralist,xcolor,xspace,amssymb}
\usepackage{times}
\newcommand{\swift}{{\sc Swift}\xspace}
%opening
\title{FMM in SWIFT}
\author{Matthieu Schaller}
\begin{document}
\maketitle
\input{gravity_derivatives}
\bibliographystyle{mnras}
\bibliography{./bibliography.bib}
\end{document}
We use the multi-index notation of \cite{Dehnen2014} to simplify expressions.
\subsection{Derivatives of the gravitational potential}
The calculation of all the
$D_\mathbf{n}(x,y,z) \equiv \nabla^{\mathbf{n}}\phi(x,y,z)$ terms up
to the relevent order can be quite tedious and it is beneficial to
automatize the whole setup. Ideally, one would like to have an
expression for each of this term that is only made of multiplications
and additions of each of the coordinates and the inverse distance. We
achieve this by writing $\phi$ as a composition of functions
$\phi(u(x,y,z))$ and apply the \textit{Fa\`a di Bruno}
formula \citep[i.e. the ``chain rule'' for higher order derivatives,
e.g.][]{Hardy2006} to construct our terms:
\begin{equation}
\label{eq:faa_di_bruno}
\frac{\partial^n}{\partial x_1 \cdots \partial x_n} \phi(u)
= \sum_{A} \phi^{(|A|)}(u) \prod_{B \in
A} \frac{\partial^{|B|}}{\prod_{c\in B}\partial x_c} u(x,y,z),
\end{equation}
where $A$ is the set of all partitions of $\lbrace1,\cdots, n\rbrace$,
$B$ is a block of a partition $A$ and $|\cdot|$ denotes the
cardinality of a set. For generic functions $\phi$ and $u$ this
formula yields an untracktable number of terms; an 8th-order
derivative will have $4140$ (!) terms in the sum\footnote{The number
of terms in the sum is given by the Bell number of the same order}. \\
We choose to write
\begin{align}
\phi(x,y,z) &= 1 / \sqrt{u(x,y,z)}, \\
u(x,y,z) &= x^2 + y^2 + z^2.
\end{align}
This choice allows to have derivatives of any order of $\phi(u)$ that
only depend on powers of $u$:
\begin{equation}
f^{(n)}(u) = \frac{\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2} -
n)}\frac{1}{u^{n+\frac{1}{2}}}.
\end{equation}
More importantly, this choice of decomposition allows us to have
derivatives of $u$ only up to second order in $x$, $y$ or $z$. The
number of non-zero terms in eq. \ref{eq:faa_di_bruno} is hence
drastically reduced. For instance, when computing
$D_{(4,1,3)} \equiv \frac{\partial^8}{\partial x^4 \partial y \partial
z^3} \phi$, $4100$ of the $4140$ terms will involve at least one
zero-valued derivative (e.g. $\partial^3/\partial x^3$ or
$\partial^2/\partial x\partial y$) of $u$. Furthermore, among the 40
remaining terms, many will involve the same derivatives and can be
grouped together, leaving us with a sum of six products of $x$,$y$ and
$z$. This is generally the case for most of the $D_\mathbf{n}$'s and
figuring out which terms are identical in a given set of partitions of
$\lbrace1,\cdots, n\rbrace$ is an interesting exercise in
combinatorics left for the reader \citep[see also][]{Hardy2006}.
#!/bin/bash
python kernels.py
pdflatex -jobname=fmm fmm_standalone.tex
bibtex fmm.aux
pdflatex -jobname=fmm fmm_standalone.tex
pdflatex -jobname=fmm fmm_standalone.tex
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment