Commit 2e7b35b1 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Don't initialise multipoles if they don't exist.

parent b3d0dc61
......@@ -289,31 +289,30 @@ 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 multipole_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;
}
#if 0
......
......@@ -496,7 +496,8 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (!cell_is_active(c, e)) return;
/* Reset the gravity acceleration tensors */
if (c->multipole) multipole_init(c->multipole);
if (e->policy & engine_policy_self_gravity)
multipole_init(c->multipole);
/* Recurse? */
if (c->split) {
......
......@@ -53,6 +53,42 @@ void runner_do_grav_up(struct runner *r, struct cell *c) {
/* } */
}
/**
* @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;
multipole_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.
......@@ -541,7 +577,7 @@ static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
// 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);
}
}
......
......@@ -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,
......
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