diff --git a/configure.ac b/configure.ac index 84548743a97423946c38b99e4811afff74bac45a..c33b52e83f724508c42a4ea1a1ce2235e1cbd90e 100644 --- a/configure.ac +++ b/configure.ac @@ -602,10 +602,10 @@ if test "$enable_warn" != "no"; then # We will do this by hand instead and only default to the macro for unknown compilers case "$ax_cv_c_compiler_vendor" in gnu | clang) - CFLAGS="$CFLAGS -Wall -Wextra -Wno-unused-parameter" + CFLAGS="$CFLAGS -Wall -Wextra -Wno-unused-parameter -Wshadow" ;; intel) - CFLAGS="$CFLAGS -w2 -Wunused-variable" + CFLAGS="$CFLAGS -w2 -Wunused-variable -Wshadow" ;; *) AX_CFLAGS_WARN_ALL diff --git a/src/cell.c b/src/cell.c index 4b344d475482549c1168a32b5740b86d3a8cfad4..7061996abe0f6417746fdcbc4d818f0755141c42 100644 --- a/src/cell.c +++ b/src/cell.c @@ -1917,46 +1917,28 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { /* If the foreign cell is active, we want its ti_end values. */ if (cell_is_active(ci, e)) scheduler_activate(s, ci->recv_ti); - /* Look for the local cell cj's send tasks. */ + /* Is the foreign cell active and will need stuff from us? */ if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = cj->send_xv; l != NULL && l->t->cj->nodeID != ci->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_xv task."); - scheduler_activate(s, l->t); + + scheduler_activate_send(s, cj->send_xv, ci->nodeID); /* Drift the cell which will be sent; note that not all sent particles will be drifted, only those that are needed. */ cell_activate_drift_part(cj, s); + /* If the local cell is also active, more stuff will be needed. */ if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_rho task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, cj->send_rho, ci->nodeID); #ifdef EXTRA_HYDRO_LOOP - for (l = cj->send_gradient; - l != NULL && l->t->cj->nodeID != ci->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_gradient task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, cj->send_gradient, ci->nodeID); #endif } } /* If the local cell is active, send its ti_end values. */ - if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_ti task."); - scheduler_activate(s, l->t); - } + if (cell_is_active(cj, e)) + scheduler_activate_send(s, cj->send_ti, ci->nodeID); } else if (cj->nodeID != engine_rank) { @@ -1974,47 +1956,29 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { /* If the foreign cell is active, we want its ti_end values. */ if (cell_is_active(cj, e)) scheduler_activate(s, cj->recv_ti); - /* Look for the local cell ci's send tasks. */ + /* Is the foreign cell active and will need stuff from us? */ if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = ci->send_xv; l != NULL && l->t->cj->nodeID != cj->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_xv task."); - scheduler_activate(s, l->t); + + scheduler_activate_send(s, ci->send_xv, cj->nodeID); /* Drift the cell which will be sent; note that not all sent particles will be drifted, only those that are needed. */ cell_activate_drift_part(ci, s); + /* If the local cell is also active, more stuff will be needed. */ if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_rho task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, ci->send_rho, cj->nodeID); #ifdef EXTRA_HYDRO_LOOP - for (l = ci->send_gradient; - l != NULL && l->t->cj->nodeID != cj->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_gradient task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, ci->send_gradient, cj->nodeID); #endif } } /* If the local cell is active, send its ti_end values. */ - if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_ti task."); - scheduler_activate(s, l->t); - } + if (cell_is_active(ci, e)) + scheduler_activate_send(s, ci->send_ti, cj->nodeID); } #endif } diff --git a/src/engine.c b/src/engine.c index 93c430d611cf573c643b4cf94325fb97333381a7..05822f54543125b359d33d231042f9a00f2e15e1 100644 --- a/src/engine.c +++ b/src/engine.c @@ -708,11 +708,11 @@ void engine_redistribute(struct engine *e) { size_t total = 0, g_total = 0, s_total = 0; size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0; for (int p = 0, r = 0; p < nr_nodes; p++) { - for (int s = 0; s < nr_nodes; s++) { + for (int n = 0; n < nr_nodes; n++) { total += counts[r]; g_total += g_counts[r]; s_total += s_counts[r]; - if (p == s) { + if (p == n) { unmoved += counts[r]; g_unmoved += g_counts[r]; s_unmoved += s_counts[r]; @@ -2754,47 +2754,30 @@ void engine_marktasks_mapper(void *map_data, int num_elements, /* If the foreign cell is active, we want its ti_end values. */ if (cell_is_active(ci, e)) scheduler_activate(s, ci->recv_ti); - /* Look for the local cell cj's send tasks. */ + /* Is the foreign cell active and will need stuff from us? */ if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = cj->send_xv; l != NULL && l->t->cj->nodeID != ci->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_xv task."); - scheduler_activate(s, l->t); + + struct link *l = + scheduler_activate_send(s, cj->send_xv, ci->nodeID); /* Drift the cell which will be sent at the level at which it is sent, i.e. drift the cell specified in the send task (l->t) itself. */ cell_activate_drift_part(l->t->ci, s); + /* If the local cell is also active, more stuff will be needed. */ if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = cj->send_rho; - l != NULL && l->t->cj->nodeID != ci->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_rho task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, cj->send_rho, ci->nodeID); #ifdef EXTRA_HYDRO_LOOP - for (l = cj->send_gradient; - l != NULL && l->t->cj->nodeID != ci->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_gradient task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, cj->send_gradient, ci->nodeID); #endif } } /* If the local cell is active, send its ti_end values. */ - if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = cj->send_ti; l != NULL && l->t->cj->nodeID != ci->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_ti task."); - scheduler_activate(s, l->t); - } + if (cell_is_active(cj, e)) + scheduler_activate_send(s, cj->send_ti, ci->nodeID); } else if (cj->nodeID != engine_rank) { @@ -2812,48 +2795,31 @@ void engine_marktasks_mapper(void *map_data, int num_elements, /* If the foreign cell is active, we want its ti_end values. */ if (cell_is_active(cj, e)) scheduler_activate(s, cj->recv_ti); - /* Look for the local cell ci's send tasks. */ + /* Is the foreign cell active and will need stuff from us? */ if (cell_is_active(cj, e)) { - struct link *l = NULL; - for (l = ci->send_xv; l != NULL && l->t->cj->nodeID != cj->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_xv task."); - scheduler_activate(s, l->t); + + struct link *l = + scheduler_activate_send(s, ci->send_xv, cj->nodeID); /* Drift the cell which will be sent at the level at which it is sent, i.e. drift the cell specified in the send task (l->t) itself. */ cell_activate_drift_part(l->t->ci, s); + /* If the local cell is also active, more stuff will be needed. */ if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = ci->send_rho; - l != NULL && l->t->cj->nodeID != cj->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_rho task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, ci->send_rho, cj->nodeID); #ifdef EXTRA_HYDRO_LOOP - for (l = ci->send_gradient; - l != NULL && l->t->cj->nodeID != cj->nodeID; l = l->next) - ; - if (l == NULL) error("Missing link to send_gradient task."); - scheduler_activate(s, l->t); + scheduler_activate_send(s, ci->send_gradient, cj->nodeID); #endif } } /* If the local cell is active, send its ti_end values. */ - if (cell_is_active(ci, e)) { - struct link *l = NULL; - for (l = ci->send_ti; l != NULL && l->t->cj->nodeID != cj->nodeID; - l = l->next) - ; - if (l == NULL) error("Missing link to send_ti task."); - scheduler_activate(s, l->t); - } + if (cell_is_active(ci, e)) + scheduler_activate_send(s, ci->send_ti, cj->nodeID); } #endif } diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h index 811d6fc8f902530840bcce4cf378c72ce25d0f4f..1f4310eea49a1b1b298fe00f00e1ae0802d9b688 100644 --- a/src/gravity/Default/gravity_iact.h +++ b/src/gravity/Default/gravity_iact.h @@ -95,7 +95,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated( } else { - const float r = r2 * r_inv; const float ui = r * h_inv; float W_ij; diff --git a/src/hydro/Shadowswift/voronoi3d_algorithm.h b/src/hydro/Shadowswift/voronoi3d_algorithm.h index 13242ad167f1936d12714af918bce7f41ac77335..353ef6f6e24356048d5f665fe1bcf81c4e6d6a50 100644 --- a/src/hydro/Shadowswift/voronoi3d_algorithm.h +++ b/src/hydro/Shadowswift/voronoi3d_algorithm.h @@ -1612,83 +1612,87 @@ __attribute__((always_inline)) INLINE void voronoi_intersect( continue; } if (c->orders[v] == 2) { - int j = voronoi_get_edge(c, v, 0); - int k = voronoi_get_edge(c, v, 1); - int b = voronoi_get_edgeindex(c, v, 1); - int l = 0; - safewhile(l < c->orders[j] && voronoi_get_edge(c, j, l) != k) { ++l; } - if (l == c->orders[j]) { + int jj = voronoi_get_edge(c, v, 0); + int kk = voronoi_get_edge(c, v, 1); + int bb = voronoi_get_edgeindex(c, v, 1); + int ll = 0; + safewhile(ll < c->orders[jj] && voronoi_get_edge(c, jj, ll) != kk) { + ++ll; + } + if (ll == c->orders[jj]) { int a = voronoi_get_edgeindex(c, v, 0); - /* j and k are not joined together. Replace their edges pointing to v - with a new edge pointing from j to k */ - voronoi_set_edge(c, j, a, k); - voronoi_set_edgeindex(c, j, a, b); - voronoi_set_edge(c, k, b, j); - voronoi_set_edgeindex(c, k, b, a); + /* jj and kk are not joined together. Replace their edges pointing to v + with a new edge pointing from jj to kk */ + voronoi_set_edge(c, jj, a, k); + voronoi_set_edgeindex(c, jj, a, bb); + voronoi_set_edge(c, kk, bb, jj); + voronoi_set_edgeindex(c, kk, bb, a); /* no new elements added to the stack: decrease the counter */ --low_order_index; } else { - /* just remove the edges from j to v and from k to v: create two new + /* just remove the edges from jj to v and from kk to v: create two new vertices */ - /* vertex j */ + /* vertex jj */ vindex = c->nvert; ++c->nvert; - c->vertices[3 * vindex] = c->vertices[3 * j]; - c->vertices[3 * vindex + 1] = c->vertices[3 * j + 1]; - c->vertices[3 * vindex + 2] = c->vertices[3 * j + 2]; - c->orders[vindex] = c->orders[j] - 1; + c->vertices[3 * vindex] = c->vertices[3 * jj]; + c->vertices[3 * vindex + 1] = c->vertices[3 * jj + 1]; + c->vertices[3 * vindex + 2] = c->vertices[3 * jj + 2]; + c->orders[vindex] = c->orders[jj] - 1; c->offsets[vindex] = c->offsets[vindex - 1] + c->orders[vindex - 1]; int m = 0; - for (int n = 0; n < c->orders[j]; ++n) { - int l = voronoi_get_edge(c, j, n); - if (l != v) { + for (int n = 0; n < c->orders[jj]; ++n) { + int lll = voronoi_get_edge(c, jj, n); + if (lll != v) { /* make a new edge */ - voronoi_set_edge(c, vindex, m, l); - voronoi_set_edgeindex(c, vindex, m, voronoi_get_edgeindex(c, j, n)); + voronoi_set_edge(c, vindex, m, lll); + voronoi_set_edgeindex(c, vindex, m, + voronoi_get_edgeindex(c, jj, n)); /* update the other vertex */ - voronoi_set_edge(c, l, voronoi_get_edgeindex(c, j, n), vindex); - voronoi_set_edgeindex(c, l, voronoi_get_edgeindex(c, j, n), m); + voronoi_set_edge(c, lll, voronoi_get_edgeindex(c, jj, n), vindex); + voronoi_set_edgeindex(c, lll, voronoi_get_edgeindex(c, jj, n), m); /* copy ngb information */ - voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, j, n)); + voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, jj, n)); ++m; } /* remove the old vertex */ - voronoi_set_edge(c, j, n, -1); - voronoi_set_edgeindex(c, j, n, -1); + voronoi_set_edge(c, jj, n, -1); + voronoi_set_edgeindex(c, jj, n, -1); } - /* vertex k */ + /* vertex kk */ vindex = c->nvert; ++c->nvert; - c->vertices[3 * vindex] = c->vertices[3 * k]; - c->vertices[3 * vindex + 1] = c->vertices[3 * k + 1]; - c->vertices[3 * vindex + 2] = c->vertices[3 * k + 2]; - c->orders[vindex] = c->orders[k] - 1; + c->vertices[3 * vindex] = c->vertices[3 * kk]; + c->vertices[3 * vindex + 1] = c->vertices[3 * kk + 1]; + c->vertices[3 * vindex + 2] = c->vertices[3 * kk + 2]; + c->orders[vindex] = c->orders[kk] - 1; c->offsets[vindex] = c->offsets[vindex - 1] + c->orders[vindex - 1]; m = 0; - for (int n = 0; n < c->orders[k]; ++n) { - int l = voronoi_get_edge(c, k, n); - if (l != v) { + for (int n = 0; n < c->orders[kk]; ++n) { + int lll = voronoi_get_edge(c, kk, n); + if (lll != v) { /* make a new edge */ - voronoi_set_edge(c, vindex, m, l); - voronoi_set_edgeindex(c, vindex, m, voronoi_get_edgeindex(c, k, n)); + voronoi_set_edge(c, vindex, m, lll); + voronoi_set_edgeindex(c, vindex, m, + voronoi_get_edgeindex(c, kk, n)); /* update the other vertex */ - voronoi_set_edge(c, l, voronoi_get_edgeindex(c, k, n), vindex); - voronoi_set_edgeindex(c, l, voronoi_get_edgeindex(c, k, n), m); + voronoi_set_edge(c, lll, voronoi_get_edgeindex(c, kk, n), vindex); + voronoi_set_edgeindex(c, lll, voronoi_get_edgeindex(c, kk, n), m); /* copy ngb information */ /* this one is special: we copy the ngb corresponding to the deleted edge and skip the one after that */ - if (n == b + 1) { - voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, k, b)); + if (n == bb + 1) { + voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, kk, bb)); } else { - voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, k, n)); + voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, kk, n)); } ++m; } /* remove the old vertex */ - voronoi_set_edge(c, k, n, -1); - voronoi_set_edgeindex(c, k, n, -1); + voronoi_set_edge(c, kk, n, -1); + voronoi_set_edgeindex(c, kk, n, -1); } - /* check if j or k has become an order 2 vertex */ + /* check if jj or kk has become an order 2 vertex */ /* if they have become an order 1 vertex, they were already an order 2 vertex, and they should already be in the list... */ if (c->orders[vindex] == 2) { @@ -1716,32 +1720,32 @@ __attribute__((always_inline)) INLINE void voronoi_intersect( voronoi_set_edge(c, v, 1, -1); voronoi_set_edgeindex(c, v, 1, -1); } else if (c->orders[v] == 1) { - int j = voronoi_get_edge(c, v, 0); + int jj = voronoi_get_edge(c, v, 0); /* we have to remove the edge between j and v. We create a new vertex */ vindex = c->nvert; ++c->nvert; - c->vertices[3 * vindex] = c->vertices[3 * j]; - c->vertices[3 * vindex + 1] = c->vertices[3 * j + 1]; - c->vertices[3 * vindex + 2] = c->vertices[3 * j + 2]; + c->vertices[3 * vindex] = c->vertices[3 * jj]; + c->vertices[3 * vindex + 1] = c->vertices[3 * jj + 1]; + c->vertices[3 * vindex + 2] = c->vertices[3 * jj + 2]; c->orders[vindex] = c->orders[j] - 1; c->offsets[vindex] = c->offsets[vindex - 1] + c->orders[vindex - 1]; int m = 0; - for (int k = 0; k < c->orders[j]; ++k) { - int l = voronoi_get_edge(c, j, k); - if (l != v) { + for (int kk = 0; kk < c->orders[j]; ++kk) { + int ll = voronoi_get_edge(c, jj, kk); + if (ll != v) { /* make a new edge */ - voronoi_set_edge(c, vindex, m, l); - voronoi_set_edgeindex(c, vindex, m, voronoi_get_edgeindex(c, j, k)); + voronoi_set_edge(c, vindex, m, ll); + voronoi_set_edgeindex(c, vindex, m, voronoi_get_edgeindex(c, jj, kk)); /* update the other vertex */ - voronoi_set_edge(c, l, voronoi_get_edgeindex(c, j, k), vindex); - voronoi_set_edgeindex(c, l, voronoi_get_edgeindex(c, j, k), m); + voronoi_set_edge(c, ll, voronoi_get_edgeindex(c, jj, kk), vindex); + voronoi_set_edgeindex(c, ll, voronoi_get_edgeindex(c, jj, kk), m); /* copy ngb information */ - voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, j, k)); + voronoi_set_ngb(c, vindex, m, voronoi_get_ngb(c, jj, kk)); ++m; } /* remove the old vertex */ - voronoi_set_edge(c, j, k, -1); - voronoi_set_edgeindex(c, j, k, -1); + voronoi_set_edge(c, jj, kk, -1); + voronoi_set_edgeindex(c, jj, kk, -1); } /* if the new vertex is a new order 2 vertex, add it to the stack */ if (c->orders[vindex] == 2) { diff --git a/src/minmax.h b/src/minmax.h index 9d92cd71d849dba615fdb05bc342014e0593d989..90dd87968a94d9601a87fd3b826000c166a98966 100644 --- a/src/minmax.h +++ b/src/minmax.h @@ -43,4 +43,32 @@ _a > _b ? _a : _b; \ }) +/** + * @brief Minimum of three numbers + * + * This macro evaluates its arguments exactly once. + */ +#define min3(x, y, z) \ + ({ \ + const __typeof__(x) _x = (x); \ + const __typeof__(y) _y = (y); \ + const __typeof__(z) _z = (z); \ + const __typeof__(x) _temp = min(_x, _y); \ + min(_temp, _z); \ + }) + +/** + * @brief Maximum of three numbers + * + * This macro evaluates its arguments exactly once. + */ +#define max3(x, y, z) \ + ({ \ + const __typeof__(x) _x = (x); \ + const __typeof__(y) _y = (y); \ + const __typeof__(z) _z = (z); \ + const __typeof__(x) _temp = max(_x, _y); \ + max(_temp, _z); \ + }) + #endif /* SWIFT_MINMAX_H */ diff --git a/src/multipole.h b/src/multipole.h index e408e5b6e0b38f724648e3a9bbade30b76e09db0..370b7d52d11f040fa053e45dd69ee4f2416eb0b2 100644 --- a/src/multipole.h +++ b/src/multipole.h @@ -989,11 +989,11 @@ INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga, * * Corresponds to equation (28c). * - * @param m The #multipole (content will be overwritten). + * @param multi The #multipole (content will be overwritten). * @param gparts The #gpart. * @param gcount The number of particles. */ -INLINE static void gravity_P2M(struct gravity_tensors *m, +INLINE static void gravity_P2M(struct gravity_tensors *multi, const struct gpart *gparts, int gcount) { /* Temporary variables */ @@ -1153,96 +1153,96 @@ INLINE static void gravity_P2M(struct gravity_tensors *m, #endif /* Store the data on the multipole. */ - m->m_pole.M_000 = mass; - m->r_max = sqrt(r_max2); - m->CoM[0] = com[0]; - m->CoM[1] = com[1]; - m->CoM[2] = com[2]; - m->m_pole.vel[0] = vel[0]; - m->m_pole.vel[1] = vel[1]; - m->m_pole.vel[2] = vel[2]; + multi->m_pole.M_000 = mass; + multi->r_max = sqrt(r_max2); + multi->CoM[0] = com[0]; + multi->CoM[1] = com[1]; + multi->CoM[2] = com[2]; + multi->m_pole.vel[0] = vel[0]; + multi->m_pole.vel[1] = vel[1]; + multi->m_pole.vel[2] = vel[2]; #if SELF_GRAVITY_MULTIPOLE_ORDER > 0 /* 1st order terms */ - m->m_pole.M_100 = M_100; - m->m_pole.M_010 = M_010; - m->m_pole.M_001 = M_001; + multi->m_pole.M_100 = M_100; + multi->m_pole.M_010 = M_010; + multi->m_pole.M_001 = M_001; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 1 /* 2nd order terms */ - m->m_pole.M_200 = M_200; - m->m_pole.M_020 = M_020; - m->m_pole.M_002 = M_002; - m->m_pole.M_110 = M_110; - m->m_pole.M_101 = M_101; - m->m_pole.M_011 = M_011; + multi->m_pole.M_200 = M_200; + multi->m_pole.M_020 = M_020; + multi->m_pole.M_002 = M_002; + multi->m_pole.M_110 = M_110; + multi->m_pole.M_101 = M_101; + multi->m_pole.M_011 = M_011; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 2 /* 3rd order terms */ - m->m_pole.M_300 = M_300; - m->m_pole.M_030 = M_030; - m->m_pole.M_003 = M_003; - m->m_pole.M_210 = M_210; - m->m_pole.M_201 = M_201; - m->m_pole.M_120 = M_120; - m->m_pole.M_021 = M_021; - m->m_pole.M_102 = M_102; - m->m_pole.M_012 = M_012; - m->m_pole.M_111 = M_111; + multi->m_pole.M_300 = M_300; + multi->m_pole.M_030 = M_030; + multi->m_pole.M_003 = M_003; + multi->m_pole.M_210 = M_210; + multi->m_pole.M_201 = M_201; + multi->m_pole.M_120 = M_120; + multi->m_pole.M_021 = M_021; + multi->m_pole.M_102 = M_102; + multi->m_pole.M_012 = M_012; + multi->m_pole.M_111 = M_111; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 3 /* 4th order terms */ - m->m_pole.M_400 = M_400; - m->m_pole.M_040 = M_040; - m->m_pole.M_004 = M_004; - m->m_pole.M_310 = M_310; - m->m_pole.M_301 = M_301; - m->m_pole.M_130 = M_130; - m->m_pole.M_031 = M_031; - m->m_pole.M_103 = M_103; - m->m_pole.M_013 = M_013; - m->m_pole.M_220 = M_220; - m->m_pole.M_202 = M_202; - m->m_pole.M_022 = M_022; - m->m_pole.M_211 = M_211; - m->m_pole.M_121 = M_121; - m->m_pole.M_112 = M_112; + multi->m_pole.M_400 = M_400; + multi->m_pole.M_040 = M_040; + multi->m_pole.M_004 = M_004; + multi->m_pole.M_310 = M_310; + multi->m_pole.M_301 = M_301; + multi->m_pole.M_130 = M_130; + multi->m_pole.M_031 = M_031; + multi->m_pole.M_103 = M_103; + multi->m_pole.M_013 = M_013; + multi->m_pole.M_220 = M_220; + multi->m_pole.M_202 = M_202; + multi->m_pole.M_022 = M_022; + multi->m_pole.M_211 = M_211; + multi->m_pole.M_121 = M_121; + multi->m_pole.M_112 = M_112; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 4 /* 5th order terms */ - m->m_pole.M_005 = M_005; - m->m_pole.M_014 = M_014; - m->m_pole.M_023 = M_023; - m->m_pole.M_032 = M_032; - m->m_pole.M_041 = M_041; - m->m_pole.M_050 = M_050; - m->m_pole.M_104 = M_104; - m->m_pole.M_113 = M_113; - m->m_pole.M_122 = M_122; - m->m_pole.M_131 = M_131; - m->m_pole.M_140 = M_140; - m->m_pole.M_203 = M_203; - m->m_pole.M_212 = M_212; - m->m_pole.M_221 = M_221; - m->m_pole.M_230 = M_230; - m->m_pole.M_302 = M_302; - m->m_pole.M_311 = M_311; - m->m_pole.M_320 = M_320; - m->m_pole.M_401 = M_401; - m->m_pole.M_410 = M_410; - m->m_pole.M_500 = M_500; + multi->m_pole.M_005 = M_005; + multi->m_pole.M_014 = M_014; + multi->m_pole.M_023 = M_023; + multi->m_pole.M_032 = M_032; + multi->m_pole.M_041 = M_041; + multi->m_pole.M_050 = M_050; + multi->m_pole.M_104 = M_104; + multi->m_pole.M_113 = M_113; + multi->m_pole.M_122 = M_122; + multi->m_pole.M_131 = M_131; + multi->m_pole.M_140 = M_140; + multi->m_pole.M_203 = M_203; + multi->m_pole.M_212 = M_212; + multi->m_pole.M_221 = M_221; + multi->m_pole.M_230 = M_230; + multi->m_pole.M_302 = M_302; + multi->m_pole.M_311 = M_311; + multi->m_pole.M_320 = M_320; + multi->m_pole.M_401 = M_401; + multi->m_pole.M_410 = M_410; + multi->m_pole.M_500 = M_500; #endif #if SELF_GRAVITY_MULTIPOLE_ORDER > 5 #error "Missing implementation for order >5" #endif #ifdef SWIFT_DEBUG_CHECKS - m->m_pole.num_gpart = gcount; + multi->m_pole.num_gpart = gcount; #endif } diff --git a/src/runner.c b/src/runner.c index 69c1512479e07da0aacad0a9e28bcaa6aafce104..0ceba7b04e1f76080deb14f2de3e69c41f8d174c 100644 --- a/src/runner.c +++ b/src/runner.c @@ -324,11 +324,10 @@ void runner_check_sorts(struct cell *c, int flags) { void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup, int clock) { - struct entry *finger; struct entry *fingers[8]; - struct part *parts = c->parts; - struct xpart *xparts = c->xparts; const int count = c->count; + const struct part *parts = c->parts; + struct xpart *xparts = c->xparts; float buff[8]; TIMER_TIC; @@ -427,7 +426,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup, } /* For each entry in the new sort list. */ - finger = c->sort[j]; + struct entry *finger = c->sort[j]; for (int ind = 0; ind < count; ind++) { /* Copy the minimum into the new sort array. */ @@ -506,7 +505,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int cleanup, /* Verify the sorting. */ for (int j = 0; j < 13; j++) { if (!(flags & (1 << j))) continue; - finger = c->sort[j]; + struct entry *finger = c->sort[j]; for (int k = 1; k < count; k++) { if (finger[k].d < finger[k - 1].d) error("Sorting failed, ascending array."); diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index dbf2311839f62ec25ebba95bc092d9c2306b4dea..ca814ed4bc29a1a79571627704635e2994ffc220 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -1219,17 +1219,19 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) { /* Let's check whether we need to still operate on this pair */ /* Get the distance between the CoMs at the last rebuild*/ - double dx = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0]; - double dy = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1]; - double dz = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2]; + double dx_rebuild = CoM_rebuild_i[0] - multi_j->CoM_rebuild[0]; + double dy_rebuild = CoM_rebuild_i[1] - multi_j->CoM_rebuild[1]; + double dz_rebuild = CoM_rebuild_i[2] - multi_j->CoM_rebuild[2]; /* Apply BC */ if (periodic) { - dx = nearest(dx, dim[0]); - dy = nearest(dy, dim[1]); - dz = nearest(dz, dim[2]); + dx_rebuild = nearest(dx_rebuild, dim[0]); + dy_rebuild = nearest(dy_rebuild, dim[1]); + dz_rebuild = nearest(dz_rebuild, dim[2]); } - const double r2_rebuild = dx * dx + dy * dy + dz * dz; + const double r2_rebuild = dx_rebuild * dx_rebuild + + dy_rebuild * dy_rebuild + + dz_rebuild * dz_rebuild; /* Is the criterion violated now but was OK at the last rebuild ? */ if (gravity_M2L_accept(multi_i->r_max_rebuild, multi_j->r_max_rebuild, diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c index 2f0b311218c86a60c5c0853c555cb48479249bce..421108d382c7b8108b6c44e6b7d268df041bd558 100644 --- a/src/runner_doiact_vec.c +++ b/src/runner_doiact_vec.c @@ -211,13 +211,13 @@ __attribute__((always_inline)) INLINE static void storeInteractions( vec_init_mask_true(int_mask2); /* Perform interactions. */ - for (int pjd = 0; pjd < icount_align; pjd += (NUM_VEC_PROC * VEC_SIZE)) { + for (int j = 0; j < icount_align; j += (NUM_VEC_PROC * VEC_SIZE)) { runner_iact_nonsym_2_vec_density( - &int_cache->r2q[pjd], &int_cache->dxq[pjd], &int_cache->dyq[pjd], - &int_cache->dzq[pjd], v_hi_inv, v_vix, v_viy, v_viz, - &int_cache->vxq[pjd], &int_cache->vyq[pjd], &int_cache->vzq[pjd], - &int_cache->mq[pjd], rhoSum, rho_dhSum, wcountSum, wcount_dhSum, - div_vSum, curlvxSum, curlvySum, curlvzSum, int_mask, int_mask2, 0); + &int_cache->r2q[j], &int_cache->dxq[j], &int_cache->dyq[j], + &int_cache->dzq[j], v_hi_inv, v_vix, v_viy, v_viz, &int_cache->vxq[j], + &int_cache->vyq[j], &int_cache->vzq[j], &int_cache->mq[j], rhoSum, + rho_dhSum, wcountSum, wcount_dhSum, div_vSum, curlvxSum, curlvySum, + curlvzSum, int_mask, int_mask2, 0); } /* Reset interaction count. */ diff --git a/src/scheduler.h b/src/scheduler.h index 0ef9426ed3248742de619d64f53bdd6b9954ed2c..c5ccbf43f3048dfb47d2d3eb7a7db6634b646700 100644 --- a/src/scheduler.h +++ b/src/scheduler.h @@ -109,7 +109,7 @@ struct scheduler { /* Inlined functions (for speed). */ /** - * @brief Add a task to the list of active tasks. + * @brief Add a regular task to the list of active tasks. * * @param s The #scheduler. * @param t The task to be added. @@ -123,6 +123,27 @@ __attribute__((always_inline)) INLINE static void scheduler_activate( } } +/** + * @brief Search and add an MPI send task to the list of active tasks. + * + * @param s The #scheduler. + * @param link The first element in the linked list of links for the task of + * interest. + * @param nodeID The nodeID of the foreign cell. + * + * @return The #link to the MPI send task. + */ +__attribute__((always_inline)) INLINE static struct link * +scheduler_activate_send(struct scheduler *s, struct link *link, int nodeID) { + + struct link *l = NULL; + for (l = link; l != NULL && l->t->cj->nodeID != nodeID; l = l->next) + ; + if (l == NULL) error("Missing link to send task."); + scheduler_activate(s, l->t); + return l; +} + /* Function prototypes. */ void scheduler_clear_active(struct scheduler *s); void scheduler_init(struct scheduler *s, struct space *space, int nr_tasks, diff --git a/src/space.c b/src/space.c index 52a34248cd9bf38e03e476e4937fa601b0ee9222..739955ebafc33aa05b67d0b9b5d778fccc2e0de7 100644 --- a/src/space.c +++ b/src/space.c @@ -403,7 +403,7 @@ void space_regrid(struct space *s, int verbose) { s->width[k] = s->dim[k] / cdim[k]; s->iwidth[k] = 1.0 / s->width[k]; } - const float dmin = min(s->width[0], min(s->width[1], s->width[2])); + const float dmin = min3(s->width[0], s->width[1], s->width[2]); /* Allocate the highest level of cells. */ s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2]; @@ -2693,14 +2693,14 @@ void space_init(struct space *s, const struct swift_params *params, } /* Decide on the minimal top-level cell size */ - const double dmax = max(max(s->dim[0], s->dim[1]), s->dim[2]); + const double dmax = max3(s->dim[0], s->dim[1], s->dim[2]); int maxtcells = parser_get_opt_param_int(params, "Scheduler:max_top_level_cells", space_max_top_level_cells_default); s->cell_min = 0.99 * dmax / maxtcells; /* Check that it is big enough. */ - const double dmin = min(min(s->dim[0], s->dim[1]), s->dim[2]); + const double dmin = min3(s->dim[0], s->dim[1], s->dim[2]); int needtcells = 3 * dmax / dmin; if (maxtcells < needtcells) error( diff --git a/src/statistics.c b/src/statistics.c index b06f087c70877fa5dda451300fcb4bf665b011be..4ec798aefda082b279742f643cce6e23649bc069 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -295,9 +295,9 @@ void stats_collect(const struct space *s, struct statistics *stats) { } /** - * @brief Apply final opetations on the #stats. + * @brief Apply final opetations on the #statistics. * - * @param stats The #stats to work on. + * @param stats The #statistics to work on. */ void stats_finalize(struct statistics *stats) { diff --git a/src/timestep.h b/src/timestep.h index 99ca1b4f90cc7b8894d4dfe0e0d2670da8f01d65..efc1b09aefef21add64ec1331fabd88e95fd0712 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( } /* Final time-step is minimum of hydro and gravity */ - float new_dt = min(min(new_dt_hydro, new_dt_cooling), new_dt_grav); + float new_dt = min3(new_dt_hydro, new_dt_cooling, new_dt_grav); /* Limit change in h */ const float dt_h_change = diff --git a/tests/test125cells.c b/tests/test125cells.c index 87325b4e9a2e254a2549b93b1031bf8af4538295..9e19d8889c4aab7c3ac5b02da9409901dc8a825a 100644 --- a/tests/test125cells.c +++ b/tests/test125cells.c @@ -632,8 +632,8 @@ int main(int argc, char *argv[]) { /* Reset particles. */ for (int i = 0; i < 125; ++i) { - for (int n = 0; n < cells[i]->count; ++n) - hydro_init_part(&cells[i]->parts[n], &space.hs); + for (int pid = 0; pid < cells[i]->count; ++pid) + hydro_init_part(&cells[i]->parts[pid], &space.hs); } /* First, sort stuff */ @@ -742,8 +742,8 @@ int main(int argc, char *argv[]) { } for (int i = 0; i < 125; ++i) { - for (int n = 0; n < cells[i]->count; ++n) - hydro_init_part(&cells[i]->parts[n], &space.hs); + for (int pid = 0; pid < cells[i]->count; ++pid) + hydro_init_part(&cells[i]->parts[pid], &space.hs); } } diff --git a/tests/testInteractions.c b/tests/testInteractions.c index 5bf036b17dce6755fc2b6f77aaea64851d613bb1..864cdd6aaf8c6cdaed6c4d6462cdaed3f6ab1a5b 100644 --- a/tests/testInteractions.c +++ b/tests/testInteractions.c @@ -251,7 +251,7 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, float vjzq[count] __attribute__((aligned(array_align))); /* Call serial interaction a set number of times. */ - for (int k = 0; k < runs; k++) { + for (int r = 0; r < runs; r++) { /* Reset particle to initial setup */ pi_serial = test_part; for (size_t i = 0; i < count; i++) pj_serial[i] = parts[i]; @@ -285,7 +285,7 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, dump_indv_particle_fields(serial_filename, &pj_serial[i]); /* Call vector interaction a set number of times. */ - for (int k = 0; k < runs; k++) { + for (int r = 0; r < runs; r++) { /* Reset particle to initial setup */ pi_vec = test_part; for (size_t i = 0; i < count; i++) pj_vec[i] = parts[i]; @@ -293,21 +293,22 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, /* Setup arrays for vector interaction. */ for (size_t i = 0; i < count; i++) { /* Compute the pairwise distance. */ - float r2 = 0.0f; - float dx[3]; + float my_r2 = 0.0f; + float my_dx[3]; for (int k = 0; k < 3; k++) { - dx[k] = pi_vec.x[k] - pj_vec[i].x[k]; - r2 += dx[k] * dx[k]; + my_dx[k] = pi_vec.x[k] - pj_vec[i].x[k]; + my_r2 += my_dx[k] * my_dx[k]; } - r2q[i] = r2; - dxq[i] = dx[0]; + r2q[i] = my_r2; + dxq[i] = my_dx[0]; + dyq[i] = my_dx[1]; + dzq[i] = my_dx[2]; + hiq[i] = pi_vec.h; piq[i] = &pi_vec; pjq[i] = &pj_vec[i]; - dyq[i] = dx[1]; - dzq[i] = dx[2]; mjq[i] = pj_vec[i].mass; vixq[i] = pi_vec.v[0]; viyq[i] = pi_vec.v[1]; @@ -354,17 +355,17 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, &curlvxSum, &curlvySum, &curlvzSum, mask, mask2, 0); } else { /* Only use one vector for interaction. */ - vector r2, dx, dy, dz; - r2.v = vec_load(&(r2q[i])); - dx.v = vec_load(&(dxq[i])); - dy.v = vec_load(&(dyq[i])); - dz.v = vec_load(&(dzq[i])); + vector my_r2, my_dx, my_dy, my_dz; + my_r2.v = vec_load(&(r2q[i])); + my_dx.v = vec_load(&(dxq[i])); + my_dy.v = vec_load(&(dyq[i])); + my_dz.v = vec_load(&(dzq[i])); runner_iact_nonsym_1_vec_density( - &r2, &dx, &dy, &dz, (hi_inv_vec), (vix_vec), (viy_vec), (viz_vec), - &(vjxq[i]), &(vjyq[i]), &(vjzq[i]), &(mjq[i]), &rhoSum, &rho_dhSum, - &wcountSum, &wcount_dhSum, &div_vSum, &curlvxSum, &curlvySum, - &curlvzSum, mask); + &my_r2, &my_dx, &my_dy, &my_dz, (hi_inv_vec), (vix_vec), (viy_vec), + (viz_vec), &(vjxq[i]), &(vjyq[i]), &(vjzq[i]), &(mjq[i]), &rhoSum, + &rho_dhSum, &wcountSum, &wcount_dhSum, &div_vSum, &curlvxSum, + &curlvySum, &curlvzSum, mask); } } @@ -465,13 +466,13 @@ void test_force_interactions(struct part test_part, struct part *parts, float cjq[count] __attribute__((aligned(array_align))); /* Call serial interaction a set number of times. */ - for (int k = 0; k < runs; k++) { + for (int r = 0; r < runs; r++) { /* Reset particle to initial setup */ pi_serial = test_part; for (size_t i = 0; i < count; i++) pj_serial[i] = parts[i]; /* Only dump data on first run. */ - if (k == 0) { + if (r == 0) { /* Dump state of particles before serial interaction. */ dump_indv_particle_fields(serial_filename, &pi_serial); for (size_t i = 0; i < count; i++) @@ -511,7 +512,7 @@ void test_force_interactions(struct part test_part, struct part *parts, dump_indv_particle_fields(serial_filename, &pj_serial[i]); /* Call vector interaction a set number of times. */ - for (int k = 0; k < runs; k++) { + for (int r = 0; r < runs; r++) { /* Reset particle to initial setup */ pi_vec = test_part; for (size_t i = 0; i < count; i++) pj_vec[i] = parts[i]; @@ -519,20 +520,20 @@ void test_force_interactions(struct part test_part, struct part *parts, /* Setup arrays for vector interaction. */ for (size_t i = 0; i < count; i++) { /* Compute the pairwise distance. */ - float r2 = 0.0f; - float dx[3]; + float my_r2 = 0.0f; + float my_dx[3]; for (int k = 0; k < 3; k++) { - dx[k] = pi_vec.x[k] - pj_vec[i].x[k]; - r2 += dx[k] * dx[k]; + my_dx[k] = pi_vec.x[k] - pj_vec[i].x[k]; + my_r2 += my_dx[k] * my_dx[k]; } piq[i] = &pi_vec; pjq[i] = &pj_vec[i]; - r2q[i] = r2; - dxq[i] = dx[0]; - dyq[i] = dx[1]; - dzq[i] = dx[2]; + r2q[i] = my_r2; + dxq[i] = my_dx[0]; + dyq[i] = my_dx[1]; + dzq[i] = my_dx[2]; hiq[i] = pi_vec.h; vixq[i] = pi_vec.v[0]; @@ -557,7 +558,7 @@ void test_force_interactions(struct part test_part, struct part *parts, } /* Only dump data on first run. */ - if (k == 0) { + if (r == 0) { /* Dump state of particles before vector interaction. */ dump_indv_particle_fields(vec_filename, piq[0]); for (size_t i = 0; i < count; i++) @@ -607,16 +608,16 @@ void test_force_interactions(struct part test_part, struct part *parts, &a_hydro_zSum, &h_dtSum, &v_sigSum, &entropy_dtSum, mask, mask2, 0); } else { /* Only use one vector for interaction. */ - vector r2, dx, dy, dz, hj, hj_inv; - r2.v = vec_load(&(r2q[i])); - dx.v = vec_load(&(dxq[i])); - dy.v = vec_load(&(dyq[i])); - dz.v = vec_load(&(dzq[i])); + vector my_r2, my_dx, my_dy, my_dz, hj, hj_inv; + my_r2.v = vec_load(&(r2q[i])); + my_dx.v = vec_load(&(dxq[i])); + my_dy.v = vec_load(&(dyq[i])); + my_dz.v = vec_load(&(dzq[i])); hj.v = vec_load(&hj_invq[i]); hj_inv = vec_reciprocal(hj); runner_iact_nonsym_1_vec_force( - &r2, &dx, &dy, &dz, vix_vec, viy_vec, viz_vec, rhoi_vec, + &my_r2, &my_dx, &my_dy, &my_dz, vix_vec, viy_vec, viz_vec, rhoi_vec, grad_hi_vec, pOrhoi2_vec, balsara_i_vec, ci_vec, &(vjxq[i]), &(vjyq[i]), &(vjzq[i]), &(rhojq[i]), &(grad_hjq[i]), &(pOrhoj2q[i]), &(balsarajq[i]), &(cjq[i]), &(mjq[i]), hi_inv_vec, hj_inv,