diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h index 06dbb0398287b200a59aa2bc1086e6aa4d51a420..3951f77d59e71439d2f621dd2a69fd62d39ba4f0 100644 --- a/src/gravity/Default/gravity_iact.h +++ b/src/gravity/Default/gravity_iact.h @@ -144,33 +144,37 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm( float r_x, float r_y, float r_z, float r2, float h, float h_inv, const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) { - //#if SELF_GRAVITY_MULTIPOLE_ORDER < 3 +#if SELF_GRAVITY_MULTIPOLE_ORDER < 3 float f_ij; runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000, &f_ij, pot); *f_x = -f_ij * r_x; *f_y = -f_ij * r_y; *f_z = -f_ij * r_z; -#if 0 // else +#else /* Get the inverse distance */ const float r_inv = 1.f / sqrtf(r2); - struct potential_derivatives_M2P pot; - compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &pot); + struct potential_derivatives_M2P d; + compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d); /* 1st order terms (monopole) */ - *f_x = m->M_000 * pot.D_100; - *f_y = m->M_000 * pot.D_010; - *f_z = m->M_000 * pot.D_001; + *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 * pot.D_300 + m->M_020 * pot.D_120 + m->M_002 * pot.D_102; - *f_y += m->M_200 * pot.D_210 + m->M_020 * pot.D_030 + m->M_002 * pot.D_012; - *f_z += m->M_200 * pot.D_201 + m->M_020 * pot.D_021 + m->M_002 * pot.D_003; - *f_x += m->M_110 * pot.D_210 + m->M_101 * pot.D_201 + m->M_011 * pot.D_111; - *f_y += m->M_110 * pot.D_120 + m->M_101 * pot.D_111 + m->M_011 * pot.D_021; - *f_z += m->M_110 * pot.D_111 + m->M_101 * pot.D_102 + m->M_011 * pot.D_012; + *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; + + *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; #endif } diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h index cf8aa54338b2e87e8bf5f2cc453ad7417eea5804..d698c66e418a120e5e7ebbe1cba0a5a9af8cd1f1 100644 --- a/src/gravity_derivatives.h +++ b/src/gravity_derivatives.h @@ -95,9 +95,16 @@ struct potential_derivatives_M2L { */ struct potential_derivatives_M2P { + /* 0th order terms */ + float D_000; + /* 1st order terms */ float D_100, D_010, D_001; + /* 2nd order terms */ + float D_200, D_020, D_002; + float D_110, D_101, D_011; + /* 3rd order terms */ float D_300, D_030, D_003; float D_210, D_201; @@ -368,11 +375,22 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2, const float r_y3 = r_y2 * r_y; const float r_z3 = r_z2 * r_z; + /* 0th order derivative */ + pot->D_000 = Dt_1; + /* 1st order derivatives */ pot->D_100 = r_x * Dt_3; pot->D_010 = r_y * Dt_3; pot->D_001 = r_z * Dt_3; + /* 2nd order derivatives */ + pot->D_200 = r_x2 * Dt_5 + Dt_3; + pot->D_020 = r_y2 * Dt_5 + Dt_3; + pot->D_002 = r_z2 * Dt_5 + Dt_3; + pot->D_110 = r_x * r_y * Dt_5; + pot->D_101 = r_x * r_z * Dt_5; + pot->D_011 = r_y * r_z * Dt_5; + /* 3rd order derivatives */ pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5; pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5; diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c index a463f19089d976360f90e8a980386ab031dc126b..5fbfb04ecd0b4a8c05c07d622dd77bd38af970ae 100644 --- a/tests/testPotentialPair.c +++ b/tests/testPotentialPair.c @@ -229,11 +229,11 @@ int main() { /* Reset the accelerations */ for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]); - /*********************************/ - /* Now, test the PM interactions */ - /*********************************/ + /**********************************/ + /* Test the basic PM interactions */ + /**********************************/ - /* Set an opening angle that allows M-P interactions */ + /* Set an opening angle that allows P-M interactions */ props.theta_crit2 = 1.; ci.gparts[0].mass = 0.; @@ -266,7 +266,98 @@ int main() { check_value(gp->a_grav[0], acc_true, "acceleration"); } - message("\n\t\t P-M interactions all good\n"); + message("\n\t\t basic P-M interactions all good\n"); + + /* Reset the accelerations */ + for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]); + +/***************************************/ +/* Test the high-order PM interactions */ +/***************************************/ + +#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3 + + /* Let's make ci more interesting */ + free(ci.gparts); + ci.gcount = 8; + if (posix_memalign((void **)&ci.gparts, gpart_align, + ci.gcount * sizeof(struct gpart)) != 0) + error("Error allocating gparts for cell ci"); + bzero(ci.gparts, ci.gcount * sizeof(struct gpart)); + + /* Place particles on a simple cube of side-length 0.2 */ + for (int n = 0; n < 8; ++n) { + if (n & 1) + ci.gparts[n].x[0] = 0.0 - 0.1; + else + ci.gparts[n].x[0] = 0.0 + 0.1; + + if (n & 2) + ci.gparts[n].x[1] = 0.5 - 0.1; + else + ci.gparts[n].x[1] = 0.5 + 0.1; + + if (n & 2) + ci.gparts[n].x[2] = 0.5 - 0.1; + else + ci.gparts[n].x[2] = 0.5 + 0.1; + + ci.gparts[n].mass = 1. / 8.; + + ci.gparts[n].epsilon = eps; + ci.gparts[n].time_bin = 1; + ci.gparts[n].type = swift_type_dark_matter; + ci.gparts[n].id_or_neg_offset = 1; +#ifdef SWIFT_DEBUG_CHECKS + ci.gparts[n].ti_drift = 8; +#endif + } + + /* Now let's make a multipole out of it. */ + gravity_reset(ci.multipole); + gravity_P2M(ci.multipole, ci.gparts, ci.gcount); + + // message("CoM=[%e %e %e]", ci.multipole->CoM[0], ci.multipole->CoM[1], + // ci.multipole->CoM[2]); + gravity_multipole_print(&ci.multipole->m_pole); + + /* Compute the forces */ + runner_dopair_grav_pp(&r, &ci, &cj); + + /* Verify everything */ + for (int n = 0; n < num_tests; ++n) { + const struct gpart *gp = &cj.gparts[n]; + const struct gravity_tensors *mpole = ci.multipole; + + double pot_true = 0, acc_true[3] = {0., 0., 0.}; + + for (int i = 0; i < 8; ++i) { + const struct gpart *gp2 = &ci.gparts[i]; + + const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1], + gp2->x[2] - gp->x[2]}; + const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]); + + pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX); + acc_true[0] -= + acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d; + acc_true[1] -= + acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d; + acc_true[2] -= + acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d; + } + + message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e", + gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true[0], gp->potential, + pot_true, acc_true[1], acc_true[2]); + + // check_value(gp->potential, pot_true, "potential"); + // check_value(gp->a_grav[0], acc_true[0], "acceleration"); + } + + message("\n\t\t high-order P-M interactions all good\n"); + +#endif free(ci.multipole); free(cj.multipole);