Commit 107aa7e5 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added the option to reconstruct all the multipoles every time-step (expensive)

parent 49a0a67c
......@@ -27,6 +27,7 @@ Valid options are:
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential.
-G Run with self-gravity.
-M Reconstruct the multipoles every time-step.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with hydrodynamics.
-S Run with stars.
......
......@@ -82,6 +82,8 @@ void print_help_message() {
"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", "-M", "",
"Reconstruct the multipoles every time-step.");
printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps. When unset use the time_end "
"parameter to stop.");
......@@ -164,6 +166,7 @@ int main(int argc, char *argv[]) {
int with_stars = 0;
int with_fp_exceptions = 0;
int with_drift_all = 0;
int with_mpole_reconstruction = 0;
int verbose = 0;
int nr_threads = 1;
int with_verbose_timers = 0;
......@@ -172,7 +175,8 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:Tv:y:")) != -1) switch (c) {
while ((c = getopt(argc, argv, "acCdDef:FgGhMn:sSt:Tv:y:")) != -1)
switch (c) {
case 'a':
with_aff = 1;
break;
......@@ -210,6 +214,9 @@ int main(int argc, char *argv[]) {
case 'h':
if (myrank == 0) print_help_message();
return 0;
case 'M':
with_mpole_reconstruction = 1;
break;
case 'n':
if (sscanf(optarg, "%d", &nsteps) != 1) {
if (myrank == 0) printf("Error parsing fixed number of steps.\n");
......@@ -521,6 +528,8 @@ int main(int argc, char *argv[]) {
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_drift_all) engine_policies |= engine_policy_drift_all;
if (with_mpole_reconstruction)
engine_policies |= engine_policy_reconstruct_mpoles;
if (with_hydro) engine_policies |= engine_policy_hydro;
if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
......
......@@ -1102,6 +1102,95 @@ void cell_reset_task_counters(struct cell *c) {
#endif
}
/**
* @brief Recursively construct all the multipoles in a cell hierarchy.
*
* @param c The #cell.
*/
void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
/* Reset everything */
gravity_reset(c->multipole);
if (c->split) {
/* Compute CoM of all progenies */
double CoM[3] = {0., 0., 0.};
double mass = 0.;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct gravity_tensors *m = c->progeny[k]->multipole;
CoM[0] += m->CoM[0] * m->m_pole.M_000;
CoM[1] += m->CoM[1] * m->m_pole.M_000;
CoM[2] += m->CoM[2] * m->m_pole.M_000;
mass += m->m_pole.M_000;
}
}
c->multipole->CoM[0] = CoM[0] / mass;
c->multipole->CoM[1] = CoM[1] / mass;
c->multipole->CoM[2] = CoM[2] / mass;
/* Now shift progeny multipoles and add them up */
struct multipole temp;
double r_max = 0.;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct cell *cp = c->progeny[k];
const struct multipole *m = &cp->multipole->m_pole;
/* Contribution to multipole */
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM);
gravity_multipole_add(&c->multipole->m_pole, &temp);
/* Upper limit of max CoM<->gpart distance */
const double dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
const double dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
const double dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
const double r2 = dx * dx + dy * dy + dz * dz;
r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
}
}
/* Alternative upper limit of max CoM<->gpart distance */
const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
? c->multipole->CoM[0] - c->loc[0]
: c->loc[0] + c->width[0] - c->multipole->CoM[0];
const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
? c->multipole->CoM[1] - c->loc[1]
: c->loc[1] + c->width[1] - c->multipole->CoM[1];
const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
? c->multipole->CoM[2] - c->loc[2]
: c->loc[2] + c->width[2] - c->multipole->CoM[2];
/* Take minimum of both limits */
c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
} else {
if (c->gcount > 0) {
gravity_P2M(c->multipole, c->gparts, c->gcount);
const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
? c->multipole->CoM[0] - c->loc[0]
: c->loc[0] + c->width[0] - c->multipole->CoM[0];
const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
? c->multipole->CoM[1] - c->loc[1]
: c->loc[1] + c->width[1] - c->multipole->CoM[1];
const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
? c->multipole->CoM[2] - c->loc[2]
: c->loc[2] + c->width[2] - c->multipole->CoM[2];
c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
} 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.;
c->multipole->r_max = 0.;
}
}
c->ti_old_multipole = ti_current;
}
/**
* @brief Computes the multi-pole brutally and compare to the
* recursively computed one.
......
......@@ -351,6 +351,7 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts);
int cell_link_sparts(struct cell *c, struct spart *sparts);
void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
void cell_check_particle_drift_point(struct cell *c, void *data);
......
......@@ -3225,6 +3225,15 @@ void engine_step(struct engine *e) {
/* Are we drifting everything (a la Gadget/GIZMO) ? */
if (e->policy & engine_policy_drift_all) engine_drift_all(e);
/* Are we reconstructing the multipoles or drifting them ?*/
if (e->policy & engine_policy_self_gravity) {
if (e->policy & engine_policy_reconstruct_mpoles)
engine_reconstruct_multipoles(e);
else
engine_drift_top_multipoles(e);
}
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
......@@ -3248,9 +3257,6 @@ 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);
......@@ -3447,6 +3453,39 @@ void engine_drift_top_multipoles(struct engine *e) {
clocks_getunit());
}
void engine_do_reconstruct_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) {
/* Construct the multipoles in this cell hierarchy */
cell_make_multipoles(c, e->ti_current);
}
}
}
/**
* @brief Reconstruct all the multipoles at all the levels in the tree.
*
* @param e The #engine.
*/
void engine_reconstruct_multipoles(struct engine *e) {
const ticks tic = getticks();
threadpool_map(&e->threadpool, engine_do_reconstruct_multipoles_mapper,
e->s->cells_top, e->s->nr_cells, sizeof(struct cell), 10, e);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
/**
* @brief Create and fill the proxies.
*
......
......@@ -66,9 +66,10 @@ enum engine_policy {
engine_policy_external_gravity = (1 << 9),
engine_policy_cosmology = (1 << 10),
engine_policy_drift_all = (1 << 11),
engine_policy_cooling = (1 << 12),
engine_policy_sourceterms = (1 << 13),
engine_policy_stars = (1 << 14)
engine_policy_reconstruct_mpoles = (1 << 12),
engine_policy_cooling = (1 << 13),
engine_policy_sourceterms = (1 << 14),
engine_policy_stars = (1 << 15)
};
extern const char *engine_policy_names[];
......@@ -256,6 +257,7 @@ 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_reconstruct_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,
......
......@@ -1232,12 +1232,10 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
* @param m_b The #multipole 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_M2M(struct multipole *m_a,
const struct multipole *m_b,
const double pos_a[3], const double pos_b[3],
int periodic) {
const double pos_a[3], const double pos_b[3]) {
/* Shift bulk velocity */
m_a->vel[0] = m_b->vel[0];
m_a->vel[1] = m_b->vel[1];
......
......@@ -2091,8 +2091,7 @@ void space_split_recursive(struct space *s, struct cell *c,
const struct multipole *m = &cp->multipole->m_pole;
/* Contribution to multipole */
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
s->periodic);
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM);
gravity_multipole_add(&c->multipole->m_pole, &temp);
/* Upper limit of max CoM<->gpart distance */
......
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