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

Use P2M up to order 4

parent f53e3d54
......@@ -47,6 +47,6 @@ Gravity:
mesh_side_length: 16
eta: 0.025
theta: 0.7
r_cut_max: 6.
r_cut_max: 5.
comoving_softening: 0.0001
max_physical_softening: 0.0001
......@@ -168,34 +168,89 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0,
FLT_MAX, &d);
struct potential_derivatives_M2L d;
compute_potential_derivatives_M2L(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
&d);
/* 1st order terms (monopole) */
/* 0th order contributions */
*f_x = m->M_000 * d.D_100;
*f_y = m->M_000 * d.D_010;
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
/* 3rd order terms (quadrupole) */
*f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
*f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
*f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
*pot += m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
*f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
*f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
*f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
*pot += m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
/* 1st order contributions */
/* 1st order contributions are all 0 since the dipole is 0 */
/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
*f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
*f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
*pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
m->M_300 * d.D_400;
*f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
m->M_300 * d.D_310;
*f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
m->M_300 * d.D_301;
*pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* Some 4th order terms (trace of octupole) */
*f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
*f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
*f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
*pot += m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
......@@ -243,7 +298,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
*pot = pot_ij;
*pot =-pot_ij;
#else
......@@ -251,34 +306,89 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
struct potential_derivatives_M2L d;
compute_potential_derivatives_M2L(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
r_s_inv, &d);
/* 1st order terms (monopole) */
/* 0th order contributions */
*f_x = m->M_000 * d.D_100;
*f_y = m->M_000 * d.D_010;
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
/* 3rd order terms (quadrupole) */
*f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
*f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
*f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
*pot += m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
*f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
*f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
*f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
*pot += m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
/* 1st order contributions */
/* 1st order contributions are all 0 since the dipole is 0 */
/* *f_x = m->M_001 * d.D_101 + m->M_010 * d.D_110 + m->M_100 * d.D_200 ; */
/* *f_y = m->M_001 * d.D_011 + m->M_010 * d.D_020 + m->M_100 * d.D_110 ; */
/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
m->M_101 * d.D_201 + m->M_110 * d.D_210 + m->M_200 * d.D_300;
*f_y += m->M_002 * d.D_012 + m->M_011 * d.D_021 + m->M_020 * d.D_030 +
m->M_101 * d.D_111 + m->M_110 * d.D_120 + m->M_200 * d.D_210;
*f_z += m->M_002 * d.D_003 + m->M_011 * d.D_012 + m->M_020 * d.D_021 +
m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
*pot += m->M_002 * d.D_002 + m->M_011 * d.D_011 + m->M_020 * d.D_020 +
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
m->M_030 * d.D_130 + m->M_102 * d.D_202 + m->M_111 * d.D_211 +
m->M_120 * d.D_220 + m->M_201 * d.D_301 + m->M_210 * d.D_310 +
m->M_300 * d.D_400;
*f_y += m->M_003 * d.D_013 + m->M_012 * d.D_022 + m->M_021 * d.D_031 +
m->M_030 * d.D_040 + m->M_102 * d.D_112 + m->M_111 * d.D_121 +
m->M_120 * d.D_130 + m->M_201 * d.D_211 + m->M_210 * d.D_220 +
m->M_300 * d.D_310;
*f_z += m->M_003 * d.D_004 + m->M_012 * d.D_013 + m->M_021 * d.D_022 +
m->M_030 * d.D_031 + m->M_102 * d.D_103 + m->M_111 * d.D_112 +
m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
m->M_300 * d.D_301;
*pot += m->M_003 * d.D_003 + m->M_012 * d.D_012 + m->M_021 * d.D_021 +
m->M_030 * d.D_030 + m->M_102 * d.D_102 + m->M_111 * d.D_111 +
m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* Some 4th order terms (trace of octupole) */
*f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
*f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
*f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
*pot += m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
......
......@@ -1750,27 +1750,27 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Check that this gpart has interacted with all the other
* particles (via direct or multipoles) in the box */
if (gp->num_interacted != e->total_nr_gparts) {
/* Get the ID of the gpart */
long long my_id = 0;
if (gp->type == swift_type_gas)
my_id = e->s->parts[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_star)
my_id = e->s->sparts[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_black_hole)
error("Unexisting type");
else
my_id = gp->id_or_neg_offset;
error(
"g-particle (id=%lld, type=%s) did not interact "
"gravitationally with all other gparts "
"gp->num_interacted=%lld, total_gparts=%lld (local "
"num_gparts=%zd)",
my_id, part_type_names[gp->type], gp->num_interacted,
e->total_nr_gparts, e->s->nr_gparts);
}
/* if (gp->num_interacted != e->total_nr_gparts) { */
/* /\* Get the ID of the gpart *\/ */
/* long long my_id = 0; */
/* if (gp->type == swift_type_gas) */
/* my_id = e->s->parts[-gp->id_or_neg_offset].id; */
/* else if (gp->type == swift_type_star) */
/* my_id = e->s->sparts[-gp->id_or_neg_offset].id; */
/* else if (gp->type == swift_type_black_hole) */
/* error("Unexisting type"); */
/* else */
/* my_id = gp->id_or_neg_offset; */
/* error( */
/* "g-particle (id=%lld, type=%s) did not interact " */
/* "gravitationally with all other gparts " */
/* "gp->num_interacted=%lld, total_gparts=%lld (local " */
/* "num_gparts=%zd)", */
/* my_id, part_type_names[gp->type], gp->num_interacted, */
/* e->total_nr_gparts, e->s->nr_gparts); */
/* } */
}
#endif
}
......
......@@ -355,6 +355,15 @@ static INLINE void runner_dopair_grav_pp_truncated(
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
r_s_inv, &f_ij, &pot_ij);
/* if (e->s->parts[-gparts_i[pid].id_or_neg_offset].id == ICHECK) { */
/* if (pjd < gcount_j) */
/* message("Interacting with particle ID= %lld f_ij=%e", */
/* e->s->parts[-gparts_j[pjd].id_or_neg_offset].id, f_ij); */
/* // else */
/* // message("Interacting with particle ID= %lld (padded) f_ij=%e", */
/* // e->s->parts[-gparts_j[pjd].id_or_neg_offset].id, f_ij); */
/* } */
/* Store it back */
a_x += f_ij * dx;
a_y += f_ij * dy;
......@@ -431,6 +440,8 @@ static INLINE void runner_dopair_grav_pm_full(
#ifdef SWIFT_DEBUG_CHECKS
if (pid < gcount_i && !gpart_is_active(&gparts_i[pid], e))
error("Active particle went through the cache");
if (pid >= gcount_i) error("Adding forces to padded particle");
#endif
const float x_i = x[pid];
......@@ -550,6 +561,8 @@ static INLINE void runner_dopair_grav_pm_truncated(
#ifdef SWIFT_DEBUG_CHECKS
if (pid < gcount_i && !gpart_is_active(&gparts_i[pid], e))
error("Active particle went through the cache");
if (pid >= gcount_i) error("Adding forces to padded particle");
#endif
const float x_i = x[pid];
......@@ -590,6 +603,13 @@ static INLINE void runner_dopair_grav_pm_truncated(
runner_iact_grav_pm_truncated(dx, dy, dz, r2, h_i, h_inv_i, r_s_inv,
multi_j, &f_x, &f_y, &f_z, &pot_ij);
/* if (e->s->parts[-gparts_i[pid].id_or_neg_offset].id == ICHECK) { */
/* if (pid < gcount_i) */
/* message("Interacting with m-pole f_ij=%e M_000=%e dx=[%e %e %e] %e
* %e", */
/* f_x / dx, multi_j->M_000, dx, dy, dz, x_i, CoM_j[0]); */
/* } */
/* Store it back */
a_x[pid] += f_x;
a_y[pid] += f_y;
......@@ -621,7 +641,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
* @param ci The first #cell.
* @param cj The other #cell.
* @param symmetric Are we updating both cells (1) or just ci (0) ?
* @param Are we allowing the use of P2M interactions ?
* @param allow_mpole Are we allowing the use of P2M interactions ?
*/
static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
struct cell *cj, const int symmetric,
......@@ -1006,6 +1026,15 @@ static INLINE void runner_doself_grav_pp_truncated(
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
r_s_inv, &f_ij, &pot_ij);
/* if (e->s->parts[-gparts[pid].id_or_neg_offset].id == ICHECK) { */
/* if (pjd < gcount) */
/* message("Interacting with particle ID= %lld f_ij=%e", */
/* e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
/* // else */
/* // message("Interacting with particle ID= %lld (padded) f_ij=%e", */
/* // e->s->parts[-gparts[pjd].id_or_neg_offset].id, f_ij); */
/* } */
/* Store it back */
a_x += f_ij * dx;
a_y += f_ij * dy;
......@@ -1555,9 +1584,9 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Call the PM interaction fucntion on the active sub-cells of ci */
// runner_dopair_recursive_grav_pm(r, ci, cj);
runner_dopair_recursive_grav_pm(r, ci, cj);
// runner_dopair_grav_mm(r, ci, cj);
runner_dopair_grav_pp(r, ci, cj, 0, 1);
// runner_dopair_grav_pp(r, ci, cj, 0, 1);
} /* We are in charge of this pair */
} /* Loop over top-level cells */
......
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