Commit 6f774986 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Implemented potential calculation for M2P interaction up to quadrupole order....

Implemented potential calculation for M2P interaction up to quadrupole order. Still need to check accuracy.
parent e14c543e
......@@ -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
}
......
......@@ -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;
......
......@@ -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);
......
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