Commit 8a46038c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Zero the gravity caches before using them. Add to the caches instead of overwritting them.

parent 7e6c42de
......@@ -46,6 +46,7 @@ Cosmology:
Gravity:
mesh_side_length: 16
eta: 0.025
theta: 0.8
theta: 0.5
r_cut_max: 6.
comoving_softening: 0.0001
max_physical_softening: 0.0001
......@@ -149,6 +149,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
const float h, const float h_inv, const struct multipole *m, float *f_x,
float *f_y, float *f_z, float *pot) {
error("www");
/* In the case where the order is < 3, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
......@@ -235,14 +237,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* In the case where the order is < 3, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if 1 // SELF_GRAVITY_MULTIPOLE_ORDER < 3
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
m->M_000, r_s_inv, &f_ij, &pot_ij);
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
*f_x = -f_ij * r_x;
*f_y = -f_ij * r_y;
*f_z = -f_ij * r_z;
*pot = pot_ij;
#else
......
......@@ -135,6 +135,35 @@ static INLINE void gravity_cache_init(struct gravity_cache *c,
c->count = padded_count;
}
/**
* @param Zero all the output fields (acceleration and potential) of a
* #gravity_cache.
*
* @param c The #gravity_cache to zero.
* @param gcount_padded The padded size of the cache arrays.
*/
__attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
struct gravity_cache *c, const int gcount_padded) {
#ifdef SWIFT_DEBUG_CHECKS
if (gcount_padded % VEC_SIZE != 0)
error("Padded gcount size not a multiple of the vector length");
#endif
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, a_x, c->a_x, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_y, c->a_y, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, a_z, c->a_z, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, pot, c->pot, SWIFT_CACHE_ALIGNMENT);
swift_assume_size(gcount_padded, VEC_SIZE);
/* Zero everything */
bzero(a_x, gcount_padded * sizeof(float));
bzero(a_y, gcount_padded * sizeof(float));
bzero(a_z, gcount_padded * sizeof(float));
bzero(pot, gcount_padded * sizeof(float));
}
/**
* @brief Fills a #gravity_cache structure with some #gpart and shift them.
*
......@@ -222,6 +251,9 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
active[i] = 0;
use_mpole[i] = 0;
}
/* Zero the output as well */
gravity_cache_zero_output(c, gcount_padded);
}
/**
......@@ -284,6 +316,9 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
m[i] = 0.f;
active[i] = 0;
}
/* Zero the output as well */
gravity_cache_zero_output(c, gcount_padded);
}
/**
......@@ -301,12 +336,16 @@ gravity_cache_populate_no_mpole(const timebin_t max_active_bin,
*/
__attribute__((always_inline)) INLINE static void
gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
const int periodic, const float dim[3],
struct gravity_cache *c,
const struct gpart *restrict gparts,
const int gcount, const int gcount_padded,
const struct cell *cell,
const struct cell *cell, const float CoM[3],
const float r_max2,
const struct gravity_props *grav_props) {
const float theta_crit2 = grav_props->theta_crit2;
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
swift_declare_aligned_ptr(float, y, c->y, SWIFT_CACHE_ALIGNMENT);
......@@ -327,6 +366,24 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
m[i] = gparts[i].mass;
active[i] = (int)(gparts[i].time_bin <= max_active_bin);
use_mpole[i] = 1;
#ifdef SWIFT_DEBUG_CHECKS
/* Distance to the CoM of the other cell. */
float dx = x[i] - CoM[0];
float dy = y[i] - CoM[1];
float dz = z[i] - CoM[2];
/* Apply periodic BC */
if (periodic) {
dx = nearestf(dx, dim[0]);
dy = nearestf(dy, dim[1]);
dz = nearestf(dz, dim[2]);
}
const float r2 = dx * dx + dy * dy + dz * dz;
if (!gravity_M2P_accept(r_max2, theta_crit2, r2))
error("Using m-pole where the test fails");
#endif
}
#ifdef SWIFT_DEBUG_CHECKS
......@@ -350,6 +407,9 @@ gravity_cache_populate_all_mpole(const timebin_t max_active_bin,
active[i] = 0;
use_mpole[i] = 0;
}
/* Zero the output as well */
gravity_cache_zero_output(c, gcount_padded);
}
/**
......
......@@ -2427,7 +2427,7 @@ __attribute__((always_inline)) INLINE static int gravity_M2P_accept(
// MATTHIEU: Make this mass-dependent ?
/* Multipole acceptance criterion (Dehnen 2002, eq.10) */
return (r2 * theta_crit2 * 0.1 > r_max2);
return (r2 * theta_crit2 > r_max2);
}
#endif /* SWIFT_MULTIPOLE_H */
......@@ -92,7 +92,7 @@ static INLINE void runner_do_grav_down(struct runner *r, struct cell *c,
/* Leaf case */
/* We can abort early if no interactions via multipole happened */
if (!c->multipole->pot.interacted) return;
// if (!c->multipole->pot.interacted) return;
if (!cell_are_gpart_drifted(c, e)) error("Un-drifted gparts");
......@@ -242,10 +242,10 @@ static INLINE void runner_dopair_grav_pp_full(
}
/* Store everything back in cache */
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
ci_cache->a_x[pid] += a_x;
ci_cache->a_y[pid] += a_y;
ci_cache->a_z[pid] += a_z;
ci_cache->pot[pid] += pot;
}
}
......@@ -368,10 +368,10 @@ static INLINE void runner_dopair_grav_pp_truncated(
}
/* Store everything back in cache */
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
ci_cache->a_x[pid] += a_x;
ci_cache->a_y[pid] += a_y;
ci_cache->a_z[pid] += a_z;
ci_cache->pot[pid] += pot;
}
}
......@@ -474,10 +474,10 @@ static INLINE void runner_dopair_grav_pm_full(
&f_z, &pot_ij);
/* Store it back */
a_x[pid] = f_x;
a_y[pid] = f_y;
a_z[pid] = f_z;
pot[pid] = pot_ij;
a_x[pid] += f_x;
a_y[pid] += f_y;
a_z[pid] += f_z;
pot[pid] += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter */
......@@ -572,7 +572,7 @@ static INLINE void runner_dopair_grav_pm_truncated(
const float r2 = dx * dx + dy * dy + dz * dz;
#ifdef SWIFT_DEBUG_CHECKSa
#ifdef SWIFT_DEBUG_CHECKS
const float r_max_j = cj->multipole->r_max;
const float r_max2 = r_max_j * r_max_j;
const float theta_crit2 = e->gravity_properties->theta_crit2;
......@@ -591,10 +591,10 @@ static INLINE void runner_dopair_grav_pm_truncated(
multi_j, &f_x, &f_y, &f_z, &pot_ij);
/* Store it back */
a_x[pid] = f_x;
a_y[pid] = f_y;
a_z[pid] = f_z;
pot[pid] = pot_ij;
a_x[pid] += f_x;
a_y[pid] += f_y;
a_z[pid] += f_z;
pot[pid] += pot_ij;
#ifdef SWIFT_DEBUG_CHECKS
/* Update the interaction counter */
......@@ -692,11 +692,11 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
/* Fill the caches */
gravity_cache_populate(e->max_active_bin, periodic, dim, ci_cache, ci->gparts,
gcount_i, gcount_padded_i, shift_i, CoM_j, rmax2_j, ci,
e->gravity_properties);
gcount_i, gcount_padded_i, shift_i, CoM_j,
rmax2_j, ci, e->gravity_properties);
gravity_cache_populate(e->max_active_bin, periodic, dim, cj_cache, cj->gparts,
gcount_j, gcount_padded_j, shift_j, CoM_i, rmax2_i, cj,
e->gravity_properties);
gcount_j, gcount_padded_j, shift_j, CoM_i,
rmax2_i, cj, e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
......@@ -751,9 +751,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
multi_j, dim, r_s_inv, e, ci->gparts,
gcount_i, cj);
if (0)
runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
multi_j, dim, r_s_inv, e, ci->gparts,
gcount_i, cj);
}
if (cj_active && symmetric) {
......@@ -763,9 +764,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
multi_i, dim, r_s_inv, e, cj->gparts,
gcount_j, ci);
if (0)
runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
multi_i, dim, r_s_inv, e, cj->gparts,
gcount_j, ci);
}
} else {
......@@ -781,8 +783,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
periodic, dim, e, ci->gparts, gcount_i, cj);
if (0)
runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
periodic, dim, e, ci->gparts, gcount_i,
cj);
}
if (cj_active && symmetric) {
......@@ -792,8 +796,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
periodic, dim, e, cj->gparts, gcount_j, ci);
if (0)
runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
periodic, dim, e, cj->gparts, gcount_j,
ci);
}
}
}
......@@ -899,10 +905,10 @@ static INLINE void runner_doself_grav_pp_full(
}
/* Store everything back in cache */
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
ci_cache->a_x[pid] += a_x;
ci_cache->a_y[pid] += a_y;
ci_cache->a_z[pid] += a_z;
ci_cache->pot[pid] += pot;
}
}
......@@ -1009,10 +1015,10 @@ static INLINE void runner_doself_grav_pp_truncated(
}
/* Store everything back in cache */
ci_cache->a_x[pid] = a_x;
ci_cache->a_y[pid] = a_y;
ci_cache->a_z[pid] = a_z;
ci_cache->pot[pid] = pot;
ci_cache->a_x[pid] += a_x;
ci_cache->a_y[pid] += a_y;
ci_cache->a_z[pid] += a_z;
ci_cache->pot[pid] += pot;
}
}
......@@ -1209,17 +1215,18 @@ static INLINE void runner_dopair_recursive_grav_pm(struct runner *r,
error("Not enough space in the cache! gcount_i=%d", gcount_i);
#endif
/* Fill the cache */
gravity_cache_populate_all_mpole(e->max_active_bin, ci_cache, ci->gparts,
gcount_i, gcount_padded_i, ci,
e->gravity_properties);
/* Recover the multipole info and the CoM locations */
const struct multipole *multi_j = &cj->multipole->m_pole;
const float r_max = cj->multipole->r_max;
const float CoM_j[3] = {(float)(cj->multipole->CoM[0]),
(float)(cj->multipole->CoM[1]),
(float)(cj->multipole->CoM[2])};
/* Fill the cache */
gravity_cache_populate_all_mpole(
e->max_active_bin, periodic, dim, ci_cache, ci->gparts, gcount_i,
gcount_padded_i, ci, CoM_j, r_max * r_max, e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
......@@ -1540,9 +1547,13 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
continue;
}
/* Call the PM interaction fucntion on the active sub-cells of ci */
++cc;
runner_dopair_recursive_grav_pm(r, ci, cj);
/* Call the PM interaction fucntion on the active sub-cells of ci */
// runner_dopair_recursive_grav_pm(r, ci, cj);
runner_dopair_grav_mm(r, ci, cj);
// runner_dopair_grav_pp(r, ci, cj, 0);
} /* 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