From 98196af3b637b95dfc75b35c13196cea25dbcb42 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Thu, 2 Apr 2020 11:56:53 +0200
Subject: [PATCH] In space_parts_get_cell_index_mapper() and friends, add an
additional check that a particle has not been wrapped back onto the exact
edge of the box due to rounding differences between the region around 0 and
the region around dim[:]
---
src/space.c | 52 ++++++++++++++++++++++++++++++++++++++++------------
1 file changed, 40 insertions(+), 12 deletions(-)
diff --git a/src/space.c b/src/space.c
index 9bfa984b82..2ed86a9ff1 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2205,9 +2205,16 @@ void space_parts_get_cell_index_mapper(void *map_data, int nr_parts,
#endif
/* Put it back into the simulation volume */
- const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
- const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
- const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+ double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Treat the case where a particle was wrapped back exactly onto
+ * the edge because of rounding issues (more accuracy around 0
+ * than around dim) */
+ if (pos_x == dim_x) pos_x = 0.0;
+ if (pos_y == dim_y) pos_y = 0.0;
+ if (pos_z == dim_z) pos_z = 0.0;
/* Get its cell index */
const int index =
@@ -2338,9 +2345,16 @@ void space_gparts_get_cell_index_mapper(void *map_data, int nr_gparts,
#endif
/* Put it back into the simulation volume */
- const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
- const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
- const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+ double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Treat the case where a particle was wrapped back exactly onto
+ * the edge because of rounding issues (more accuracy around 0
+ * than around dim) */
+ if (pos_x == dim_x) pos_x = 0.0;
+ if (pos_y == dim_y) pos_y = 0.0;
+ if (pos_z == dim_z) pos_z = 0.0;
/* Get its cell index */
const int index =
@@ -2477,9 +2491,16 @@ void space_sparts_get_cell_index_mapper(void *map_data, int nr_sparts,
#endif
/* Put it back into the simulation volume */
- const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
- const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
- const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+ double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Treat the case where a particle was wrapped back exactly onto
+ * the edge because of rounding issues (more accuracy around 0
+ * than around dim) */
+ if (pos_x == dim_x) pos_x = 0.0;
+ if (pos_y == dim_y) pos_y = 0.0;
+ if (pos_z == dim_z) pos_z = 0.0;
/* Get its cell index */
const int index =
@@ -2612,9 +2633,16 @@ void space_bparts_get_cell_index_mapper(void *map_data, int nr_bparts,
#endif
/* Put it back into the simulation volume */
- const double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
- const double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
- const double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+ double pos_x = box_wrap(old_pos_x, 0.0, dim_x);
+ double pos_y = box_wrap(old_pos_y, 0.0, dim_y);
+ double pos_z = box_wrap(old_pos_z, 0.0, dim_z);
+
+ /* Treat the case where a particle was wrapped back exactly onto
+ * the edge because of rounding issues (more accuracy around 0
+ * than around dim) */
+ if (pos_x == dim_x) pos_x = 0.0;
+ if (pos_y == dim_y) pos_y = 0.0;
+ if (pos_z == dim_z) pos_z = 0.0;
/* Get its cell index */
const int index =
--
GitLab