Skip to content
Snippets Groups Projects
Commit a3fdc554 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Prepare for merge

parent 5736f328
No related branches found
No related tags found
1 merge request!324Gravity multi dt
Hide whitespace changes
Inline Side-by-side
src/multipole.h
+ 6
− 18
View file @ a3fdc554
......@@ -906,55 +906,45 @@ INLINE static void gravity_M2M(struct multipole *m_a,
m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
/* Shift 2nd order term */
m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
m_a->M_020 = m_b->M_020 + X_010(dx) * m_b->M_010 + X_020(dx) * m_b->M_000;
m_a->M_002 = m_b->M_002 + X_001(dx) * m_b->M_001 + X_002(dx) * m_b->M_000;
m_a->M_110 = m_b->M_110 + X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +
X_110(dx) * m_b->M_000;
m_a->M_101 = m_b->M_101 + X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +
X_101(dx) * m_b->M_000;
m_a->M_011 = m_b->M_011 + X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +
X_011(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* Shift 3rd order term */
m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + X_200(dx) * m_b->M_100 +
X_300(dx) * m_b->M_000;
m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + X_020(dx) * m_b->M_010 +
X_030(dx) * m_b->M_000;
m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + X_002(dx) * m_b->M_001 +
X_003(dx) * m_b->M_000;
m_a->M_210 = m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 +
X_210(dx) * m_b->M_000;
m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 +
X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +
X_201(dx) * m_b->M_000;
m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 +
X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +
X_120(dx) * m_b->M_000;
m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 +
X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +
X_021(dx) * m_b->M_000;
m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 +
X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +
X_102(dx) * m_b->M_000;
m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 +
X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +
X_012(dx) * m_b->M_000;
m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 +
X_001(dx) * m_b->M_110 + X_110(dx) * m_b->M_001 +
X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +
......@@ -1013,7 +1003,6 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b->F_010 += m_a->M_000 * D_010(dx, dy, dz, r_inv);
l_b->F_001 += m_a->M_000 * D_001(dx, dy, dz, r_inv);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order multipole term (addition to rank 0)*/
......@@ -1043,7 +1032,6 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
l_b->F_101 += m_a->M_000 * D_101(dx, dy, dz, r_inv);
l_b->F_011 += m_a->M_000 * D_011(dx, dy, dz, r_inv);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order multipole term (addition to rank 0)*/
......@@ -1156,7 +1144,6 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
la->F_010 += X_000(dx) * lb->F_010;
la->F_001 += X_000(dx) * lb->F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* Shift 2nd order multipole term (addition to rank 0)*/
......@@ -1181,7 +1168,6 @@ INLINE static void gravity_L2L(struct grav_tensor *la,
la->F_101 += X_000(dx) * lb->F_101;
la->F_011 += X_000(dx) * lb->F_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* Shift 3rd order multipole term (addition to rank 0)*/
......@@ -1271,6 +1257,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
gp->a_grav[2] += X_000(dx) * lb->F_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 1st order interaction */
gp->a_grav[0] +=
X_100(dx) * lb->F_200 + X_010(dx) * lb->F_110 + X_001(dx) * lb->F_101;
......@@ -1280,6 +1267,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
X_100(dx) * lb->F_101 + X_010(dx) * lb->F_011 + X_001(dx) * lb->F_002;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 2nd order interaction */
gp->a_grav[0] +=
X_200(dx) * lb->F_300 + X_020(dx) * lb->F_120 + X_002(dx) * lb->F_102;
......
src/runner_doiact_grav.h
+ 13
− 62
View file @ a3fdc554
......@@ -25,6 +25,14 @@
#include "gravity.h"
#include "part.h"
/**
* @brief Recursively propagate the multipoles down the tree by applying the
* L2L and L2P kernels.
*
* @param r The #runner.
* @param c The #cell we are working on.
* @param timer Are we timing this ?
*/
void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
......@@ -114,62 +122,6 @@ void runner_dopair_grav_pm(const struct runner *r,
const struct cell *restrict cj) {
error("Function should not be called");
/* const struct engine *e = r->e; */
/* const int gcount = ci->gcount; */
/* struct gpart *restrict gparts = ci->gparts; */
/* const struct gravity_tensors *multi = cj->multipole; */
/* 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 (gcount == 0) error("Empty cell!"); */
/* if (multi->m_pole.mass == 0.0) */
/* error("Multipole does not seem to have been set."); */
/* #endif */
/* /\* Anything to do here? *\/ */
/* if (!cell_is_active(ci, e)) return; */
/* #if ICHECK > 0 */
/* for (int pid = 0; pid < gcount; pid++) { */
/* /\* Get a hold of the ith part in ci. *\/ */
/* struct gpart *restrict gp = &gparts[pid]; */
/* if (gp->id_or_neg_offset == ICHECK) */
/* message("id=%lld loc=[ %f %f %f ] size= %f count= %d", */
/* gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2], */
/* cj->width[0], cj->gcount); */
/* } */
/* #endif */
/* /\* Loop over every particle in leaf. *\/ */
/* for (int pid = 0; pid < gcount; pid++) { */
/* /\* Get a hold of the ith part in ci. *\/ */
/* struct gpart *restrict gp = &gparts[pid]; */
/* if (!gpart_is_active(gp, e)) continue; */
/* /\* Compute the pairwise distance. *\/ */
/* const float dx[3] = {multi->CoM[0] - gp->x[0], // x */
/* multi->CoM[1] - gp->x[1], // y */
/* multi->CoM[2] - gp->x[2]}; // z */
/* const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
/* /\* Interact !*\/ */
/* runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole); */
/* #ifdef SWIFT_DEBUG_CHECKS */
/* gp->mass_interacted += multi->m_pole.mass; */
/* #endif */
/* } */
/* TIMER_TOC(timer_dopair_grav_pm); */
}
/**
......@@ -307,8 +259,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
if (c->gcount == 0) // MATTHIEU sanity check
error("Empty cell !");
if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
/* Anything to do here? */
......@@ -403,7 +354,8 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
const int gcount_j = cj->gcount;
/* Early abort? */
if (gcount_i == 0 || gcount_j == 0) error("Empty cell !");
if (gcount_i == 0 || gcount_j == 0)
error("Doing pair gravity on an empty cell !");
/* Bad stuff will happen if cell sizes are different */
if (ci->width[0] != cj->width[0])
......@@ -500,7 +452,7 @@ void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
#ifdef SWIFT_DEBUG_CHECKS
/* Early abort? */
if (c->gcount == 0) error("Empty cell !");
if (c->gcount == 0) error("Doing self gravity on an empty cell !");
#endif
TIMER_TIC;
......@@ -594,8 +546,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* cj->loc[1] - pos_i[1], // y */
/* cj->loc[2] - pos_i[2]}; // z */
/* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
// if (r2 > max_d2) continue;
/* if (r2 > max_d2) continue; */
if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment