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

Cleaned-up a bit the code. Calculation of the best-fitting correction to the...

Cleaned-up a bit the code. Calculation of the best-fitting correction to the bias can now be implemented.
parent 72b9580f
Branches
No related tags found
No related merge requests found
......@@ -135,7 +135,6 @@ for orientation in range( 26 ):
# continue
cell_edge = sqrt(sum(abs(array(this_axis))))# - sum(abs(this_axis[this_axis<0]))
print this_axis
# print 2*cell_edge
print "actual:", orientation, min(dist), max(dist), min_dist[0], max_dist[0]
print "---"
continue
......
......@@ -1184,7 +1184,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
struct part_local *parts_i, *parts_j;
float xi[3];
float dx[3], ai[3], mi, mj, r2, w, ir;
float min_dist, max_dist;
float min_dist, max_dist, length_dist, mid_dist;
#ifdef SANITY_CHECKS
......@@ -1194,10 +1194,14 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
#endif
/* Get the sorted indices and stuff. */
/* Get the sorted indices and geometry. */
struct index *ind_i, *ind_j;
struct multipole multi;
get_axis(&ci, &cj, &ind_i, &ind_j, &min_dist, &max_dist);
length_dist = max_dist - min_dist;
mid_dist = min_dist + 0.5f * length_dist;
/* Initialise the multipole */
struct multipole multi;
multipole_reset(&multi);
/* Allocate and fill-in the local parts. */
......@@ -1210,6 +1214,9 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
(struct part_local *)alloca(sizeof(struct part_local) * count_j)) ==
NULL)
error("Failed to allocate local part arrays.");
/* Store only the particle position relative to a point close to the cells */
/* This point does not have to be in the middle cj->loc is good enough */
float midpoint[3] = { cj->loc[0] , cj->loc[1], cj->loc[2] };
for (i = 0; i < count_i; i++) {
int pid = ind_i[i].ind;
......@@ -1251,7 +1258,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Loop over all particles in ci... */
for (i = count_i - 1; i >= 0; i--) {
/* Get a local copy of the distance along the axis. */
/* Get a local copy of the maximal distance along the axis. */
float di = parts_i[i].d + d_max;
/* Init the ith particle's data. */
......@@ -1303,7 +1310,7 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
} /* loop over every other particle. */
/* Add any remaining particles to the COM. */
/* Add any remaining particles to the multipole. */
for (int jj = j; jj < count_j; jj++) {
multipole_add(&multi, parts_j[jj].x, parts_j[jj].mass);
}
......@@ -1311,12 +1318,12 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Shrink count_j to the latest valid particle. */
count_j = j;
/* Interact part_i with the center of mass. */
/* Interact part_i with the multipole. */
multipole_iact(&multi, xi, ai);
#if ICHECK >= 0
if (parts_i[i].id == ICHECK)
message("Interaction with multipole"); //, parts_j[j].id );
message("Interaction with multipole");
#endif
#ifdef COUNTERS
......@@ -1353,12 +1360,24 @@ static inline void iact_pair_direct_sorted_multipole(struct cell *ci,
/* Interact part_j with the COM. */
multipole_iact(&multi, parts_j[j].x, parts_j[j].a);
#if ICHECK >= 0
if (parts_j[j].id == ICHECK)
message("Interaction with multipole");
#endif
#ifdef COUNTERS
++count_direct_sorted_pm_j;
#endif
}
/* Apply the correction to cancel the bias */
for (i = 0; i < count_i; i++) {
message("Length = %f Mid = %f dist = %f", length_dist, mid_dist, parts_i[i].d);
}
/* Copy the accelerations back to the original particles. */
for (i = 0; i < count_i; i++) {
int pid = ind_i[i].ind;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment