Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • dc-oman1/swiftsim
  • swift/swiftsim
  • pdraper/swiftsim
  • tkchan/swiftsim
  • dc-turn5/swiftsim
5 results
Select Git revision
Show changes
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
* Matthieu Schaller (schaller@strw.leidenuniv.nl)
* 2015 Peter W. Draper (p.w.draper@durham.ac.uk)
* Angus Lepper (angus.lepper@ed.ac.uk)
* 2016 John A. Regan (john.a.regan@durham.ac.uk)
* Tom Theuns (tom.theuns@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
/* Config parameters. */
#include <config.h>
/* Some standard headers. */
#include <stdlib.h>
#include <unistd.h>
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
#endif
/* Load the profiler header, if needed. */
#ifdef WITH_PROFILER
#include <gperftools/profiler.h>
#endif
/* This object's header. */
#include "engine.h"
/* Local headers. */
#include "active.h"
#include "atomic.h"
#include "cell.h"
#include "clocks.h"
#include "cycle.h"
#include "debug.h"
#include "error.h"
#include "feedback.h"
#include "proxy.h"
#include "timers.h"
/**
* @brief Mark tasks to be un-skipped and set the sort flags accordingly.
* Threadpool mapper function.
*
* @param map_data pointer to the tasks
* @param num_elements number of tasks
* @param extra_data pointer to int that will define if a rebuild is needed.
*/
void engine_marktasks_mapper(void *map_data, int num_elements,
void *extra_data) {
/* Unpack the arguments. */
struct task *tasks = (struct task *)map_data;
size_t *rebuild_space = &((size_t *)extra_data)[1];
struct scheduler *s = (struct scheduler *)(((size_t *)extra_data)[2]);
struct engine *e = (struct engine *)((size_t *)extra_data)[0];
const int nodeID = e->nodeID;
const int with_timestep_limiter = e->policy & engine_policy_timestep_limiter;
const int with_timestep_sync = e->policy & engine_policy_timestep_sync;
const int with_sinks = e->policy & engine_policy_sinks;
const int with_stars = e->policy & engine_policy_stars;
const int with_star_formation = (e->policy & engine_policy_star_formation);
const int with_star_formation_sink = with_sinks && with_stars;
const int with_feedback = e->policy & engine_policy_feedback;
for (int ind = 0; ind < num_elements; ind++) {
/* Get basic task information */
struct task *t = &tasks[ind];
const enum task_types t_type = t->type;
const enum task_subtypes t_subtype = t->subtype;
/* Single-cell task? */
if (t_type == task_type_self || t_type == task_type_sub_self) {
/* Local pointer. */
struct cell *ci = t->ci;
#ifdef SWIFT_DEBUG_CHECKS
if (ci->nodeID != nodeID) error("Non-local self task found");
#endif
const int ci_active_hydro = cell_is_active_hydro(ci, e);
const int ci_active_gravity = cell_is_active_gravity(ci, e);
const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
const int ci_active_sinks =
cell_is_active_sinks(ci, e) || ci_active_hydro;
const int ci_active_stars = cell_need_activating_stars(
ci, e, with_star_formation, with_star_formation_sink);
const int ci_active_rt = cell_is_rt_active(ci, e);
/* Activate the hydro drift */
if (t_type == task_type_self && t_subtype == task_subtype_density) {
if (ci_active_hydro) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
if (with_timestep_limiter) cell_activate_limiter(ci, s);
}
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_density) {
if (ci_active_hydro) {
scheduler_activate(s, t);
cell_activate_subcell_hydro_tasks(ci, NULL, s, with_timestep_limiter);
if (with_timestep_limiter) cell_activate_limiter(ci, s);
}
}
else if (t_type == task_type_self && t_subtype == task_subtype_force) {
if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_force) {
if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t->type == task_type_self &&
t->subtype == task_subtype_limiter) {
if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t->type == task_type_sub_self &&
t->subtype == task_subtype_limiter) {
if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_self && t_subtype == task_subtype_gradient) {
if (ci_active_hydro) scheduler_activate(s, t);
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_gradient) {
if (ci_active_hydro) scheduler_activate(s, t);
}
/* Activate the star density */
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_density) {
if (ci_active_stars) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
cell_activate_drift_spart(ci, s);
if (with_timestep_sync) cell_activate_sync_part(ci, s);
}
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_density) {
if (ci_active_stars) {
scheduler_activate(s, t);
cell_activate_subcell_stars_tasks(ci, NULL, s, with_star_formation,
with_star_formation_sink,
with_timestep_sync);
}
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_prep1) {
if (ci_active_stars) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_prep1) {
if (ci_active_stars) scheduler_activate(s, t);
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_prep2) {
if (ci_active_stars) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_prep2) {
if (ci_active_stars) scheduler_activate(s, t);
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_stars_feedback) {
if (ci_active_stars) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_stars_feedback) {
if (ci_active_stars) scheduler_activate(s, t);
}
/* Activate the sink formation */
else if (t_type == task_type_self &&
t_subtype == task_subtype_sink_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
cell_activate_drift_sink(ci, s);
cell_activate_sink_formation_tasks(ci->top, s);
if (with_timestep_sync) cell_activate_sync_part(ci, s);
}
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_sink_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
cell_activate_subcell_sinks_tasks(ci, NULL, s, with_timestep_sync);
}
}
/* Activate the sink merger */
else if (t_type == task_type_self &&
t_subtype == task_subtype_sink_do_sink_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_sink_do_sink_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
}
}
/* Activate the sink accretion */
else if (t_type == task_type_self &&
t_subtype == task_subtype_sink_do_gas_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_sink_do_gas_swallow) {
if (ci_active_sinks) {
scheduler_activate(s, t);
}
}
/* Activate the black hole density */
else if (t_type == task_type_self &&
t_subtype == task_subtype_bh_density) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
cell_activate_drift_part(ci, s);
cell_activate_drift_bpart(ci, s);
if (with_timestep_sync) cell_activate_sync_part(ci, s);
}
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_bh_density) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
cell_activate_subcell_black_holes_tasks(ci, NULL, s,
with_timestep_sync);
}
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_bh_swallow) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_bh_swallow) {
if (ci_active_black_holes) scheduler_activate(s, t);
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_do_gas_swallow) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_do_gas_swallow) {
if (ci_active_black_holes) scheduler_activate(s, t);
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_do_bh_swallow) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_do_bh_swallow) {
if (ci_active_black_holes) scheduler_activate(s, t);
}
else if (t_type == task_type_self &&
t_subtype == task_subtype_bh_feedback) {
if (ci_active_black_holes) {
scheduler_activate(s, t);
}
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_bh_feedback) {
if (ci_active_black_holes) scheduler_activate(s, t);
}
/* Activate the gravity drift */
else if (t_type == task_type_self && t_subtype == task_subtype_grav) {
if (ci_active_gravity) {
scheduler_activate(s, t);
cell_activate_subcell_grav_tasks(t->ci, NULL, s);
}
}
/* Activate the gravity drift */
else if (t_type == task_type_self &&
t_subtype == task_subtype_external_grav) {
if (ci_active_gravity) {
scheduler_activate(s, t);
cell_activate_drift_gpart(t->ci, s);
}
}
/* Activate RT tasks */
else if (t_type == task_type_self &&
t_subtype == task_subtype_rt_gradient) {
if (ci_active_rt) scheduler_activate(s, t);
}
else if (t_type == task_type_sub_self &&
t_subtype == task_subtype_rt_gradient) {
if (ci_active_rt) {
scheduler_activate(s, t);
cell_activate_subcell_rt_tasks(ci, NULL, s, /*sub_cycle=*/0);
}
}
else if (t_subtype == task_subtype_rt_transport) {
if (ci_active_rt) scheduler_activate(s, t);
}
#ifdef SWIFT_DEBUG_CHECKS
else {
error("Invalid task type / sub-type encountered");
}
#endif
}
/* Pair? */
else if (t_type == task_type_pair || t_type == task_type_sub_pair) {
/* Local pointers. */
struct cell *ci = t->ci;
struct cell *cj = t->cj;
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = cj->nodeID;
#else
const int ci_nodeID = nodeID;
const int cj_nodeID = nodeID;
#endif
const int ci_active_hydro = cell_is_active_hydro(ci, e);
const int cj_active_hydro = cell_is_active_hydro(cj, e);
const int ci_active_gravity = cell_is_active_gravity(ci, e);
const int cj_active_gravity = cell_is_active_gravity(cj, e);
const int ci_active_black_holes = cell_is_active_black_holes(ci, e);
const int cj_active_black_holes = cell_is_active_black_holes(cj, e);
const int ci_active_sinks =
cell_is_active_sinks(ci, e) || ci_active_hydro;
const int cj_active_sinks =
cell_is_active_sinks(cj, e) || cj_active_hydro;
const int ci_active_stars = cell_need_activating_stars(
ci, e, with_star_formation, with_star_formation_sink);
const int cj_active_stars = cell_need_activating_stars(
cj, e, with_star_formation, with_star_formation_sink);
const int ci_active_rt = cell_is_rt_active(ci, e);
const int cj_active_rt = cell_is_rt_active(cj, e);
/* Only activate tasks that involve a local active cell. */
if ((t_subtype == task_subtype_density ||
t_subtype == task_subtype_gradient ||
t_subtype == task_subtype_limiter ||
t_subtype == task_subtype_force) &&
((ci_active_hydro && ci_nodeID == nodeID) ||
(cj_active_hydro && cj_nodeID == nodeID))) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair && t_subtype == task_subtype_density) {
/* Store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Activate the hydro drift tasks. */
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
/* And the limiter */
if (ci_nodeID == nodeID && with_timestep_limiter)
cell_activate_limiter(ci, s);
if (cj_nodeID == nodeID && with_timestep_limiter)
cell_activate_limiter(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_density) {
cell_activate_subcell_hydro_tasks(t->ci, t->cj, s,
with_timestep_limiter);
}
}
/* Stars density */
else if ((t_subtype == task_subtype_stars_density) &&
(ci_active_stars || cj_active_stars) &&
(ci_nodeID == nodeID || cj_nodeID == nodeID)) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair) {
/* Add stars_in dependencies for each cell that is part of
* a pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->stars.stars_in);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->stars.stars_in);
/* Do ci */
if (ci_active_stars) {
/* stars for ci */
atomic_or(&ci->stars.requires_sorts, 1 << t->flags);
ci->stars.dx_max_sort_old = ci->stars.dx_max_sort;
/* hydro for cj */
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Activate the drift tasks. */
if (ci_nodeID == nodeID) cell_activate_drift_spart(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
if (cj_nodeID == nodeID && with_timestep_sync)
cell_activate_sync_part(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(cj, t->flags, s);
cell_activate_stars_sorts(ci, t->flags, s);
}
/* Do cj */
if (cj_active_stars) {
/* hydro for ci */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
/* stars for cj */
atomic_or(&cj->stars.requires_sorts, 1 << t->flags);
cj->stars.dx_max_sort_old = cj->stars.dx_max_sort;
/* Activate the drift tasks. */
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_spart(cj, s);
if (ci_nodeID == nodeID && with_timestep_sync)
cell_activate_sync_part(ci, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_stars_sorts(cj, t->flags, s);
}
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_stars_density) {
/* Add stars_in dependencies for each cell that is part of
* a pair/sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->stars.stars_in);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->stars.stars_in);
cell_activate_subcell_stars_tasks(ci, cj, s, with_star_formation,
with_star_formation_sink,
with_timestep_sync);
}
}
/* Stars prep1 */
else if (t_subtype == task_subtype_stars_prep1) {
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
(ci_active_stars || cj_active_stars)) {
scheduler_activate(s, t);
/* If there are active sparts in ci, activate hydro ghost in cj */
if (ci_active_stars)
scheduler_activate(s, cj->hydro.super->hydro.prep1_ghost);
/* If there are active sparts in cj, activate hydro ghost in ci */
if (cj_active_stars)
scheduler_activate(s, ci->hydro.super->hydro.prep1_ghost);
} else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
(cj_active_stars)) {
scheduler_activate(s, t);
/* If there are active sparts in cj, activate hydro ghost in ci */
scheduler_activate(s, ci->hydro.super->hydro.prep1_ghost);
} else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
(ci_active_stars)) {
scheduler_activate(s, t);
/* If there are active sparts in ci, activate hydro ghost in cj */
scheduler_activate(s, cj->hydro.super->hydro.prep1_ghost);
}
}
/* Stars prep2 */
else if (t_subtype == task_subtype_stars_prep2) {
/* We only want to activate the task if the cell is active and is
going to update some sparts on the *local* node */
if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
(ci_active_stars || cj_active_stars)) {
scheduler_activate(s, t);
} else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
(ci_active_stars)) {
scheduler_activate(s, t);
} else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
(cj_active_stars)) {
scheduler_activate(s, t);
}
}
/* Stars feedback */
else if (t_subtype == task_subtype_stars_feedback) {
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && cj_nodeID == nodeID) &&
(ci_active_stars || cj_active_stars)) {
scheduler_activate(s, t);
} else if ((ci_nodeID == nodeID && cj_nodeID != nodeID) &&
(cj_active_stars)) {
scheduler_activate(s, t);
} else if ((ci_nodeID != nodeID && cj_nodeID == nodeID) &&
(ci_active_stars)) {
scheduler_activate(s, t);
}
if (t->type == task_type_pair || t->type == task_type_sub_pair) {
/* Add stars_out dependencies for each cell that is part of
* a pair/sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->stars.stars_out);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->stars.stars_out);
}
}
/* Black_Holes density */
else if ((t_subtype == task_subtype_bh_density ||
t_subtype == task_subtype_bh_swallow ||
t_subtype == task_subtype_do_gas_swallow ||
t_subtype == task_subtype_do_bh_swallow ||
t_subtype == task_subtype_bh_feedback) &&
(ci_active_black_holes || cj_active_black_holes) &&
(ci_nodeID == nodeID || cj_nodeID == nodeID)) {
scheduler_activate(s, t);
/* Set the correct drifting flags */
if (t_type == task_type_pair && t_subtype == task_subtype_bh_density) {
if (ci_nodeID == nodeID) cell_activate_drift_bpart(ci, s);
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
if (cj_nodeID == nodeID) cell_activate_drift_bpart(cj, s);
/* Activate bh_in for each cell that is part of
* a pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->black_holes.black_holes_in);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->black_holes.black_holes_in);
}
if ((t_type == task_type_pair || t_type == task_type_sub_pair) &&
t_subtype == task_subtype_bh_feedback) {
/* Add bh_out dependencies for each cell that is part of
* a pair/sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->black_holes.black_holes_out);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->black_holes.black_holes_out);
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_bh_density) {
cell_activate_subcell_black_holes_tasks(ci, cj, s,
with_timestep_sync);
/* Activate bh_in for each cell that is part of
* a sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->black_holes.black_holes_in);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->black_holes.black_holes_in);
}
}
/* Gravity */
else if ((t_subtype == task_subtype_grav) &&
((ci_active_gravity && ci_nodeID == nodeID) ||
(cj_active_gravity && cj_nodeID == nodeID))) {
scheduler_activate(s, t);
if (t_type == task_type_pair && t_subtype == task_subtype_grav) {
/* Activate the gravity drift */
cell_activate_subcell_grav_tasks(t->ci, t->cj, s);
}
#ifdef SWIFT_DEBUG_CHECKS
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_grav) {
error("Invalid task sub-type encountered");
}
#endif
}
/* Sink formation */
else if ((t_subtype == task_subtype_sink_swallow ||
t_subtype == task_subtype_sink_do_sink_swallow ||
t_subtype == task_subtype_sink_do_gas_swallow) &&
(ci_active_sinks || cj_active_sinks) &&
(ci_nodeID == nodeID || cj_nodeID == nodeID)) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair &&
t_subtype == task_subtype_sink_swallow) {
/* Activate the sink drift for the sink merger */
if (ci_nodeID == nodeID) {
cell_activate_drift_sink(ci, s);
cell_activate_sink_formation_tasks(ci->top, s);
/* Activate all sink_in tasks for each cell involved
* in pair type tasks */
scheduler_activate(s, ci->hydro.super->sinks.sink_in);
}
if (cj_nodeID == nodeID) {
cell_activate_drift_sink(cj, s);
if (ci->top != cj->top) {
cell_activate_sink_formation_tasks(cj->top, s);
}
/* Activate all sink_in tasks for each cell involved
* in pair type tasks */
scheduler_activate(s, cj->hydro.super->sinks.sink_in);
}
/* Do ci */
if (ci_active_sinks) {
/* hydro for cj */
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Activate the drift tasks. */
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
if (cj_nodeID == nodeID && with_timestep_sync)
cell_activate_sync_part(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(cj, t->flags, s);
}
/* Do cj */
if (cj_active_sinks) {
/* hydro for ci */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
/* Activate the drift tasks. */
/* Activate the sink drift for the merger */
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (ci_nodeID == nodeID && with_timestep_sync)
cell_activate_sync_part(ci, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
}
}
else if (t_type == task_type_sub_pair &&
t_subtype == task_subtype_sink_swallow) {
/* Activate all sink_in tasks for each cell involved
* in sub_pair type tasks */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->sinks.sink_in);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->sinks.sink_in);
/* Store current values of dx_max and h_max. */
cell_activate_subcell_sinks_tasks(ci, cj, s, with_timestep_sync);
}
else if ((t_type == task_type_pair || t_type == task_type_sub_pair) &&
t_subtype == task_subtype_sink_do_gas_swallow) {
/* Activate sinks_out for each cell that is part of
* a pair/sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->sinks.sink_out);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->sinks.sink_out);
}
}
/* RT gradient and transport tasks */
else if (t_subtype == task_subtype_rt_gradient) {
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && ci_active_rt) ||
(cj_nodeID == nodeID && cj_active_rt)) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair) {
/* Store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Check the sorts and activate them if needed. */
cell_activate_rt_sorts(ci, t->flags, s);
cell_activate_rt_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair) {
cell_activate_subcell_rt_tasks(ci, cj, s, /*sub_cycle=*/0);
}
}
}
else if (t_subtype == task_subtype_rt_transport) {
/* We only want to activate the task if the cell is active and is
going to update some gas on the *local* node */
if ((ci_nodeID == nodeID && ci_active_rt) ||
(cj_nodeID == nodeID && cj_active_rt)) {
/* The gradient and transport task subtypes mirror the hydro tasks.
* Therefore all the (subcell) sorts and drifts should already have
* been activated properly in the hydro part of the activation. */
scheduler_activate(s, t);
if (t_type == task_type_pair || t_type == task_type_sub_pair) {
/* Activate transport_out for each cell that is part of
* a pair/sub_pair task as to not miss any dependencies */
if (ci_nodeID == nodeID)
scheduler_activate(s, ci->hydro.super->rt.rt_transport_out);
if (cj_nodeID == nodeID)
scheduler_activate(s, cj->hydro.super->rt.rt_transport_out);
}
}
}
/* Pair tasks between inactive local cells and active remote cells. */
if ((ci_nodeID != nodeID && cj_nodeID == nodeID && ci_active_hydro &&
!cj_active_hydro) ||
(ci_nodeID == nodeID && cj_nodeID != nodeID && !ci_active_hydro &&
cj_active_hydro)) {
#if defined(WITH_MPI) && defined(MPI_SYMMETRIC_FORCE_INTERACTION)
if (t_subtype == task_subtype_force) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair) {
/* Store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Activate the hydro drift tasks. */
if (ci_nodeID == nodeID) cell_activate_drift_part(ci, s);
if (cj_nodeID == nodeID) cell_activate_drift_part(cj, s);
/* And the limiter */
if (ci_nodeID == nodeID && with_timestep_limiter)
cell_activate_limiter(ci, s);
if (cj_nodeID == nodeID && with_timestep_limiter)
cell_activate_limiter(cj, s);
/* Check the sorts and activate them if needed. */
cell_activate_hydro_sorts(ci, t->flags, s);
cell_activate_hydro_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair) {
cell_activate_subcell_hydro_tasks(t->ci, t->cj, s,
with_timestep_limiter);
}
}
}
/* Pair tasks between inactive local cells and active remote cells. */
if ((ci_nodeID != nodeID && cj_nodeID == nodeID && ci_active_rt &&
!cj_active_rt) ||
(ci_nodeID == nodeID && cj_nodeID != nodeID && !ci_active_rt &&
cj_active_rt)) {
if (t_subtype == task_subtype_rt_transport) {
scheduler_activate(s, t);
/* Set the correct sorting flags */
if (t_type == task_type_pair) {
/* Store some values. */
atomic_or(&ci->hydro.requires_sorts, 1 << t->flags);
atomic_or(&cj->hydro.requires_sorts, 1 << t->flags);
ci->hydro.dx_max_sort_old = ci->hydro.dx_max_sort;
cj->hydro.dx_max_sort_old = cj->hydro.dx_max_sort;
/* Check the sorts and activate them if needed. */
cell_activate_rt_sorts(ci, t->flags, s);
cell_activate_rt_sorts(cj, t->flags, s);
}
/* Store current values of dx_max and h_max. */
else if (t_type == task_type_sub_pair) {
cell_activate_subcell_rt_tasks(ci, cj, s, /*sub_cycle=*/0);
}
}
#endif
}
/* Only interested in density tasks as of here. */
if (t_subtype == task_subtype_density) {
/* Too much particle movement? */
if (cell_need_rebuild_for_hydro_pair(ci, cj)) *rebuild_space = 1;
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
if (ci_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (cj_active_hydro) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
if (ci_active_hydro) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
#ifdef EXTRA_HYDRO_LOOP
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gradient);
#endif
}
}
/* If the local cell is inactive and the remote cell is active, we
* still need to receive stuff to be able to do the force interaction
* on this node as well. */
else if (ci_active_hydro) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* NOTE: (yuyttenh, 09/2022) Since the particle communications send
* over whole particles currently, just activating the gradient
* send/recieve should be enough for now. The remote active
* particles are only needed for the sorts and the flux exchange on
* the node of the inactive cell, so sending over the xv and
* gradient suffices. If at any point the commutications change, we
* should probably also send over the rho separately. */
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
#ifndef EXTRA_HYDRO_LOOP
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
#else
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gradient);
#endif
#endif
}
/* If the foreign cell is active, we want its particles for the
* limiter */
if (ci_active_hydro && with_timestep_limiter) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_limiter);
scheduler_activate_unpack(s, ci->mpi.unpack, task_subtype_limiter);
}
/* Is the foreign cell active and will need stuff from us? */
if (ci_active_hydro) {
struct link *l = scheduler_activate_send(
s, cj->mpi.send, task_subtype_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 (cj_active_hydro) {
scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
ci_nodeID);
#ifdef EXTRA_HYDRO_LOOP
scheduler_activate_send(s, cj->mpi.send, task_subtype_gradient,
ci_nodeID);
#endif
}
}
/* If the foreign cell is inactive, but the local cell is active,
* we still need to send stuff to be able to do the force interaction
* on both nodes */
else if (cj_active_hydro) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* See NOTE on line 867 */
struct link *l = scheduler_activate_send(
s, cj->mpi.send, task_subtype_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);
#ifndef EXTRA_HYDRO_LOOP
scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
ci_nodeID);
#else
scheduler_activate_send(s, cj->mpi.send, task_subtype_gradient,
ci_nodeID);
#endif
#endif
}
/* If the local cell is active, send its particles for the limiting.
*/
if (cj_active_hydro && with_timestep_limiter) {
scheduler_activate_send(s, cj->mpi.send, task_subtype_limiter,
ci_nodeID);
scheduler_activate_pack(s, cj->mpi.pack, task_subtype_limiter,
ci_nodeID);
}
/* Propagating new star counts? */
if (with_star_formation_sink) error("TODO");
if (with_star_formation && with_feedback) {
if (ci_active_hydro && ci->hydro.count > 0) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_sf_counts);
}
if (cj_active_hydro && cj->hydro.count > 0) {
scheduler_activate_send(s, cj->mpi.send, task_subtype_sf_counts,
ci_nodeID);
}
}
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active_hydro) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
if (cj_active_hydro) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
#ifdef EXTRA_HYDRO_LOOP
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gradient);
#endif
}
}
/* If the local cell is inactive and the remote cell is active, we
* still need to receive stuff to be able to do the force interaction
* on this node as well. */
else if (cj_active_hydro) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* See NOTE on line 867. */
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
#ifndef EXTRA_HYDRO_LOOP
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
#else
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gradient);
#endif
#endif
}
/* If the foreign cell is active, we want its particles for the
* limiter */
if (cj_active_hydro && with_timestep_limiter) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_limiter);
scheduler_activate_unpack(s, cj->mpi.unpack, task_subtype_limiter);
}
/* Is the foreign cell active and will need stuff from us? */
if (cj_active_hydro) {
struct link *l = scheduler_activate_send(
s, ci->mpi.send, task_subtype_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 (ci_active_hydro) {
scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
cj_nodeID);
#ifdef EXTRA_HYDRO_LOOP
scheduler_activate_send(s, ci->mpi.send, task_subtype_gradient,
cj_nodeID);
#endif
}
}
/* If the foreign cell is inactive, but the local cell is active,
* we still need to send stuff to be able to do the force interaction
* on both nodes */
else if (ci_active_hydro) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* See NOTE on line 867. */
struct link *l = scheduler_activate_send(
s, ci->mpi.send, task_subtype_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);
#ifndef EXTRA_HYDRO_LOOP
scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
cj_nodeID);
#else
scheduler_activate_send(s, ci->mpi.send, task_subtype_gradient,
cj_nodeID);
#endif
#endif
}
/* If the local cell is active, send its particles for the limiting.
*/
if (ci_active_hydro && with_timestep_limiter) {
scheduler_activate_send(s, ci->mpi.send, task_subtype_limiter,
cj_nodeID);
scheduler_activate_pack(s, ci->mpi.pack, task_subtype_limiter,
cj_nodeID);
}
/* Propagating new star counts? */
if (with_star_formation_sink) error("TODO");
if (with_star_formation && with_feedback) {
if (cj_active_hydro && cj->hydro.count > 0) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_sf_counts);
}
if (ci_active_hydro && ci->hydro.count > 0) {
scheduler_activate_send(s, ci->mpi.send, task_subtype_sf_counts,
cj_nodeID);
}
}
}
#endif
}
/* Only interested in stars_density tasks as of here. */
else if (t->subtype == task_subtype_stars_density) {
/* Too much particle movement? */
if (cell_need_rebuild_for_stars_pair(ci, cj)) *rebuild_space = 1;
if (cell_need_rebuild_for_stars_pair(cj, ci)) *rebuild_space = 1;
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
if (ci_nodeID != nodeID) {
if (cj_active_stars) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_xv);
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_part_prep1);
#endif
/* If the local cell is active, more stuff will be needed. */
scheduler_activate_send(s, cj->mpi.send, task_subtype_spart_density,
ci_nodeID);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_send(s, cj->mpi.send, task_subtype_spart_prep2,
ci_nodeID);
#endif
cell_activate_drift_spart(cj, s);
}
if (ci_active_stars) {
scheduler_activate_recv(s, ci->mpi.recv,
task_subtype_spart_density);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_spart_prep2);
#endif
/* Is the foreign cell active and will need stuff from us? */
scheduler_activate_send(s, cj->mpi.send, task_subtype_xv,
ci_nodeID);
scheduler_activate_send(s, cj->mpi.send, task_subtype_rho,
ci_nodeID);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_send(s, cj->mpi.send, task_subtype_part_prep1,
ci_nodeID);
#endif
/* 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);
}
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active_stars) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_xv);
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_part_prep1);
#endif
/* If the local cell is active, more stuff will be needed. */
scheduler_activate_send(s, ci->mpi.send, task_subtype_spart_density,
cj_nodeID);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_send(s, ci->mpi.send, task_subtype_spart_prep2,
cj_nodeID);
#endif
cell_activate_drift_spart(ci, s);
}
if (cj_active_stars) {
scheduler_activate_recv(s, cj->mpi.recv,
task_subtype_spart_density);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_spart_prep2);
#endif
/* Is the foreign cell active and will need stuff from us? */
scheduler_activate_send(s, ci->mpi.send, task_subtype_xv,
cj_nodeID);
scheduler_activate_send(s, ci->mpi.send, task_subtype_rho,
cj_nodeID);
#ifdef EXTRA_STAR_LOOPS
scheduler_activate_send(s, ci->mpi.send, task_subtype_part_prep1,
cj_nodeID);
#endif
/* 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);
}
}
#endif
}
/* Only interested in sink_swallow tasks as of here. */
else if (t->subtype == task_subtype_sink_swallow) {
/* Too much particle movement? */
if (cell_need_rebuild_for_sinks_pair(ci, cj)) *rebuild_space = 1;
if (cell_need_rebuild_for_sinks_pair(cj, ci)) *rebuild_space = 1;
#ifdef WITH_MPI
error("TODO");
#endif
}
/* Only interested in black hole density tasks as of here. */
else if (t->subtype == task_subtype_bh_density) {
/* Too much particle movement? */
if (cell_need_rebuild_for_black_holes_pair(ci, cj)) *rebuild_space = 1;
if (cell_need_rebuild_for_black_holes_pair(cj, ci)) *rebuild_space = 1;
scheduler_activate(s, ci->hydro.super->black_holes.swallow_ghost_0);
scheduler_activate(s, cj->hydro.super->black_holes.swallow_ghost_0);
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
if (ci_nodeID != nodeID) {
/* Receive the foreign parts to compute BH accretion rates and do the
* swallowing */
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rho);
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_part_swallow);
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_merger);
/* Send the local BHs to tag the particles to swallow and to do
* feedback */
scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_rho,
ci_nodeID);
scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_feedback,
ci_nodeID);
/* Drift before you send */
cell_activate_drift_bpart(cj, s);
/* Receive the foreign BHs to tag particles to swallow and for
* feedback */
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_rho);
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_bpart_feedback);
/* Send the local part information */
scheduler_activate_send(s, cj->mpi.send, task_subtype_rho, ci_nodeID);
scheduler_activate_send(s, cj->mpi.send, task_subtype_part_swallow,
ci_nodeID);
scheduler_activate_send(s, cj->mpi.send, task_subtype_bpart_merger,
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);
} else if (cj_nodeID != nodeID) {
/* Receive the foreign parts to compute BH accretion rates and do the
* swallowing */
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rho);
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_part_swallow);
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_merger);
/* Send the local BHs to tag the particles to swallow and to do
* feedback */
scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_rho,
cj_nodeID);
scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_feedback,
cj_nodeID);
/* Drift before you send */
cell_activate_drift_bpart(ci, s);
/* Receive the foreign BHs to tag particles to swallow and for
* feedback */
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_rho);
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_bpart_feedback);
/* Send the local part information */
scheduler_activate_send(s, ci->mpi.send, task_subtype_rho, cj_nodeID);
scheduler_activate_send(s, ci->mpi.send, task_subtype_part_swallow,
cj_nodeID);
scheduler_activate_send(s, ci->mpi.send, task_subtype_bpart_merger,
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);
}
#endif
}
/* Only interested in gravity tasks as of here. */
else if (t_subtype == task_subtype_grav) {
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
if (ci_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (cj_active_gravity)
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_gpart);
/* Is the foreign cell active and will need stuff from us? */
if (ci_active_gravity) {
struct link *l = scheduler_activate_send(
s, cj->mpi.send, task_subtype_gpart, 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_gpart(l->t->ci, s);
}
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active_gravity)
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_gpart);
/* Is the foreign cell active and will need stuff from us? */
if (cj_active_gravity) {
struct link *l = scheduler_activate_send(
s, ci->mpi.send, task_subtype_gpart, 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_gpart(l->t->ci, s);
}
}
#endif
} /* Only interested in RT tasks as of here. */
else if (t->subtype == task_subtype_rt_gradient) {
#ifdef WITH_MPI
/* Activate the send/recv tasks. */
if (ci_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (cj_active_rt) {
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
if (ci_active_rt) {
/* We only need updates later on if the other cell is active too
*/
scheduler_activate_recv(s, ci->mpi.recv,
task_subtype_rt_transport);
}
} else if (ci_active_rt) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* If the local cell is inactive and the remote cell is active, we
* still need to receive stuff to be able to do the force
* interaction on this node as well. */
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_gradient);
scheduler_activate_recv(s, ci->mpi.recv, task_subtype_rt_transport);
#endif
}
/* Is the foreign cell active and will need stuff from us? */
if (ci_active_rt) {
scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
ci_nodeID);
if (cj_active_rt) {
scheduler_activate_send(s, cj->mpi.send,
task_subtype_rt_transport, ci_nodeID);
}
} else if (cj_active_rt) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* If the foreign cell is inactive, but the local cell is active,
* we still need to send stuff to be able to do the force
* interaction on both nodes */
scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_gradient,
ci_nodeID);
scheduler_activate_send(s, cj->mpi.send, task_subtype_rt_transport,
ci_nodeID);
#endif
}
} else if (cj_nodeID != nodeID) {
/* If the local cell is active, receive data from the foreign cell. */
if (ci_active_rt) {
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
if (cj_active_rt) {
/* We only need updates later on if the other cell is active too
*/
scheduler_activate_recv(s, cj->mpi.recv,
task_subtype_rt_transport);
}
} else if (cj_active_rt) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* If the local cell is inactive and the remote cell is active, we
* still need to receive stuff to be able to do the force
* interaction on this node as well. */
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_gradient);
scheduler_activate_recv(s, cj->mpi.recv, task_subtype_rt_transport);
#endif
}
/* Is the foreign cell active and will need stuff from us? */
if (cj_active_rt) {
scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
cj_nodeID);
if (ci_active_rt) {
/* We only need updates later on if the other cell is active too
*/
scheduler_activate_send(s, ci->mpi.send,
task_subtype_rt_transport, cj_nodeID);
}
} else if (ci_active_rt) {
#ifdef MPI_SYMMETRIC_FORCE_INTERACTION
/* If the foreign cell is inactive, but the local cell is active,
* we still need to send stuff to be able to do the force
* interaction on both nodes */
scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_gradient,
cj_nodeID);
scheduler_activate_send(s, ci->mpi.send, task_subtype_rt_transport,
cj_nodeID);
#endif
}
}
#endif
}
}
/* End force for hydro ? */
else if (t_type == task_type_end_hydro_force) {
if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
}
/* End force for gravity ? */
else if (t_type == task_type_end_grav_force) {
if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
}
/* Activate the weighting task for neutrinos */
else if (t_type == task_type_neutrino_weight) {
if (cell_is_active_gravity(t->ci, e)) {
scheduler_activate(s, t);
}
}
/* Kick ? */
else if (t_type == task_type_kick1 || t_type == task_type_kick2) {
if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
/* Hydro ghost tasks ? */
else if (t_type == task_type_ghost || t_type == task_type_extra_ghost ||
t_type == task_type_ghost_in || t_type == task_type_ghost_out) {
if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
}
/* csds tasks ? */
else if (t->type == task_type_csds) {
if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
cell_is_active_stars(t->ci, e))
scheduler_activate(s, t);
}
/* Gravity stuff ? */
else if (t_type == task_type_grav_down ||
t_type == task_type_grav_long_range ||
t_type == task_type_init_grav ||
t_type == task_type_init_grav_out ||
t_type == task_type_drift_gpart_out ||
t_type == task_type_grav_down_in) {
if (cell_is_active_gravity(t->ci, e)) scheduler_activate(s, t);
}
/* Multipole - Multipole interaction task */
else if (t_type == task_type_grav_mm) {
/* Local pointers. */
const struct cell *ci = t->ci;
const struct cell *cj = t->cj;
#ifdef WITH_MPI
const int ci_nodeID = ci->nodeID;
const int cj_nodeID = (cj != NULL) ? cj->nodeID : -1;
#else
const int ci_nodeID = nodeID;
const int cj_nodeID = nodeID;
#endif
const int ci_active_gravity = cell_is_active_gravity_mm(ci, e);
const int cj_active_gravity = cell_is_active_gravity_mm(cj, e);
if ((ci_active_gravity && ci_nodeID == nodeID) ||
(cj_active_gravity && cj_nodeID == nodeID))
scheduler_activate(s, t);
}
/* Star ghost tasks ? */
else if (t_type == task_type_stars_ghost ||
t_type == task_type_stars_prep_ghost1 ||
t_type == task_type_hydro_prep_ghost1 ||
t_type == task_type_stars_prep_ghost2) {
if (cell_need_activating_stars(t->ci, e, with_star_formation,
with_star_formation_sink))
scheduler_activate(s, t);
}
/* Feedback implicit tasks? */
else if (t_type == task_type_stars_in || t_type == task_type_stars_out) {
if (cell_need_activating_stars(t->ci, e, with_star_formation,
with_star_formation_sink))
scheduler_activate(s, t);
}
/* Sink implicit tasks? */
else if (t_type == task_type_sink_in || t_type == task_type_sink_out ||
t_type == task_type_sink_ghost1) {
if (cell_is_active_sinks(t->ci, e) || cell_is_active_hydro(t->ci, e))
scheduler_activate(s, t);
}
/* Black hole ghost tasks ? */
else if (t_type == task_type_bh_density_ghost ||
t_type == task_type_bh_swallow_ghost1 ||
t_type == task_type_bh_swallow_ghost2 ||
t_type == task_type_bh_swallow_ghost3) {
if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
}
/* Black holes implicit tasks? */
else if (t_type == task_type_bh_in || t_type == task_type_bh_out) {
if (cell_is_active_black_holes(t->ci, e)) scheduler_activate(s, t);
}
/* Time-step or time-step collection? */
else if (t_type == task_type_timestep || t_type == task_type_collect) {
t->ci->hydro.updated = 0;
t->ci->grav.updated = 0;
t->ci->stars.updated = 0;
t->ci->sinks.updated = 0;
t->ci->black_holes.updated = 0;
if (!cell_is_empty(t->ci)) {
if (cell_is_active_hydro(t->ci, e) ||
cell_is_active_gravity(t->ci, e) ||
cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
}
else if ((t_type == task_type_send && t_subtype == task_subtype_tend) ||
(t_type == task_type_recv && t_subtype == task_subtype_tend)) {
if (!cell_is_empty(t->ci)) {
scheduler_activate(s, t);
}
}
/* Subgrid tasks: cooling */
else if (t_type == task_type_cooling || t_type == task_type_cooling_in ||
t_type == task_type_cooling_out) {
if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
}
/* Subgrid tasks: star formation */
else if (t_type == task_type_star_formation) {
if (cell_is_active_hydro(t->ci, e)) {
cell_activate_star_formation_tasks(t->ci, s, with_feedback);
cell_activate_super_spart_drifts(t->ci, s);
}
}
/* Subgrid tasks: star formation from sinks */
else if (t_type == task_type_star_formation_sink) {
if (cell_is_active_hydro(t->ci, e) || cell_is_active_sinks(t->ci, e)) {
cell_activate_star_formation_sink_tasks(t->ci, s, with_feedback);
cell_activate_super_spart_drifts(t->ci, s);
}
}
/* Radiative transfer implicit tasks */
else if (t->type == task_type_rt_in) {
if (cell_is_rt_active(t->ci, e) ||
cell_need_activating_stars(t->ci, e, with_star_formation,
with_star_formation_sink))
scheduler_activate(s, t);
}
else if (t->type == task_type_rt_ghost1 || t->type == task_type_rt_ghost2 ||
t->type == task_type_rt_transport_out ||
t->type == task_type_rt_tchem ||
t->type == task_type_rt_advance_cell_time ||
t->type == task_type_rt_out) {
if (cell_is_rt_active(t->ci, e)) scheduler_activate(s, t);
/* Note that rt_collect_times never needs to be active on main steps,
* which is always what follows engine_marktasks().*/
}
/* Subgrid tasks: sink formation */
else if (t_type == task_type_sink_formation) {
if (with_star_formation_sink && t->ci->hydro.count > 0 &&
cell_is_active_hydro(t->ci, e)) {
cell_activate_sink_formation_tasks(t->ci, s);
cell_activate_super_sink_drifts(t->ci, s);
}
}
}
}
/**
* @brief Mark tasks to be un-skipped and set the sort flags accordingly.
*
* @return 1 if the space has to be rebuilt, 0 otherwise.
*/
int engine_marktasks(struct engine *e) {
struct scheduler *s = &e->sched;
const ticks tic = getticks();
int rebuild_space = 0;
/* Run through the tasks and mark as skip or not. */
size_t extra_data[3] = {(size_t)e, (size_t)rebuild_space, (size_t)&e->sched};
threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks,
sizeof(struct task), threadpool_auto_chunk_size, extra_data);
rebuild_space = extra_data[1];
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
/* All is well... */
return rebuild_space;
}
......@@ -264,8 +264,8 @@ struct redist_mapper_data {
int *dest = \
mydata->dest + (ptrdiff_t)(parts - (struct TYPE *)mydata->base); \
int *lcounts = NULL; \
if ((lcounts = (int *)calloc( \
sizeof(int), mydata->nr_nodes * mydata->nr_nodes)) == NULL) \
if ((lcounts = (int *)calloc(mydata->nr_nodes * mydata->nr_nodes, \
sizeof(int))) == NULL) \
error("Failed to allocate counts thread-specific buffer"); \
for (int k = 0; k < num_elements; k++) { \
for (int j = 0; j < 3; j++) { \
......@@ -319,6 +319,14 @@ void ENGINE_REDISTRIBUTE_DEST_MAPPER(gpart);
*/
void ENGINE_REDISTRIBUTE_DEST_MAPPER(bpart);
/**
* @brief Accumulate the counts of sink particles per cell.
* Threadpool helper for accumulating the counts of particles per cell.
*
* sink version.
*/
void ENGINE_REDISTRIBUTE_DEST_MAPPER(sink);
#endif /* redist_mapper_data */
#ifdef WITH_MPI /* savelink_mapper_data */
......@@ -399,12 +407,22 @@ void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 1);
void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(bpart, 0);
#endif
/**
* @brief Save position of sink-gpart links.
* Threadpool helper for accumulating the counts of particles per cell.
*/
#ifdef SWIFT_DEBUG_CHECKS
void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 1);
#else
void ENGINE_REDISTRIBUTE_SAVELINK_MAPPER(sink, 0);
#endif
#endif /* savelink_mapper_data */
#ifdef WITH_MPI /* relink_mapper_data */
/* Support for relinking parts, gparts, sparts and bparts after moving between
* nodes. */
/* Support for relinking parts, gparts, sparts, bparts and sinks after moving
* between nodes. */
struct relink_mapper_data {
int nodeID;
int nr_nodes;
......@@ -412,6 +430,7 @@ struct relink_mapper_data {
int *s_counts;
int *g_counts;
int *b_counts;
int *sink_counts;
struct space *s;
};
......@@ -435,6 +454,7 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
int *g_counts = mydata->g_counts;
int *s_counts = mydata->s_counts;
int *b_counts = mydata->b_counts;
int *sink_counts = mydata->sink_counts;
struct space *s = mydata->s;
for (int i = 0; i < num_elements; i++) {
......@@ -446,12 +466,14 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
size_t offset_gparts = 0;
size_t offset_sparts = 0;
size_t offset_bparts = 0;
size_t offset_sinks = 0;
for (int n = 0; n < node; n++) {
int ind_recv = n * nr_nodes + nodeID;
offset_parts += counts[ind_recv];
offset_gparts += g_counts[ind_recv];
offset_sparts += s_counts[ind_recv];
offset_bparts += b_counts[ind_recv];
offset_sinks += sink_counts[ind_recv];
}
/* Number of gparts sent from this node. */
......@@ -493,6 +515,17 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
s->gparts[k].id_or_neg_offset = -partner_index;
s->bparts[partner_index].gpart = &s->gparts[k];
}
/* Does this gpart have a sink partner ? */
else if (s->gparts[k].type == swift_type_sink) {
const ptrdiff_t partner_index =
offset_sinks - s->gparts[k].id_or_neg_offset;
/* Re-link */
s->gparts[k].id_or_neg_offset = -partner_index;
s->sinks[partner_index].gpart = &s->gparts[k];
}
}
}
}
......@@ -519,13 +552,6 @@ void engine_redistribute_relink_mapper(void *map_data, int num_elements,
void engine_redistribute(struct engine *e) {
#ifdef WITH_MPI
#ifdef SWIFT_DEBUG_CHECKS
const int nr_sinks_new = 0;
#endif
if (e->policy & engine_policy_sinks) {
error("Not implemented yet");
}
const int nr_nodes = e->nr_nodes;
const int nodeID = e->nodeID;
struct space *s = e->s;
......@@ -536,12 +562,14 @@ void engine_redistribute(struct engine *e) {
struct gpart *gparts = s->gparts;
struct spart *sparts = s->sparts;
struct bpart *bparts = s->bparts;
struct sink *sinks = s->sinks;
ticks tic = getticks();
size_t nr_parts = s->nr_parts;
size_t nr_gparts = s->nr_gparts;
size_t nr_sparts = s->nr_sparts;
size_t nr_bparts = s->nr_bparts;
size_t nr_sinks = s->nr_sinks;
/* Start by moving inhibited particles to the end of the arrays */
for (size_t k = 0; k < nr_parts; /* void */) {
......@@ -609,6 +637,27 @@ void engine_redistribute(struct engine *e) {
}
}
/* Now move inhibited sink particles to the end of the arrays */
for (size_t k = 0; k < nr_sinks; /* void */) {
if (sinks[k].time_bin == time_bin_inhibited ||
sinks[k].time_bin == time_bin_not_created) {
nr_sinks -= 1;
/* Swap the particle */
memswap(&s->sinks[k], &s->sinks[nr_sinks], sizeof(struct sink));
/* Swap the link with the gpart */
if (s->sinks[k].gpart != NULL) {
s->sinks[k].gpart->id_or_neg_offset = -k;
}
if (s->sinks[nr_sinks].gpart != NULL) {
s->sinks[nr_sinks].gpart->id_or_neg_offset = -nr_sinks;
}
} else {
k++;
}
}
/* Finally do the same with the gravity particles */
for (size_t k = 0; k < nr_gparts; /* void */) {
if (gparts[k].time_bin == time_bin_inhibited ||
......@@ -626,6 +675,8 @@ void engine_redistribute(struct engine *e) {
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_black_hole) {
s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_sink) {
s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
if (s->gparts[nr_gparts].type == swift_type_gas) {
......@@ -637,6 +688,9 @@ void engine_redistribute(struct engine *e) {
} else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
} else if (s->gparts[nr_gparts].type == swift_type_sink) {
s->sinks[-s->gparts[nr_sinks].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
}
} else {
k++;
......@@ -648,7 +702,7 @@ void engine_redistribute(struct engine *e) {
/* Allocate temporary arrays to store the counts of particles to be sent
* and the destination of each particle */
int *counts;
if ((counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
if ((counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
error("Failed to allocate counts temporary buffer.");
int *dest;
......@@ -727,7 +781,7 @@ void engine_redistribute(struct engine *e) {
/* Get destination of each s-particle */
int *s_counts;
if ((s_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
if ((s_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
error("Failed to allocate s_counts temporary buffer.");
int *s_dest;
......@@ -793,7 +847,7 @@ void engine_redistribute(struct engine *e) {
/* Get destination of each b-particle */
int *b_counts;
if ((b_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
if ((b_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
error("Failed to allocate b_counts temporary buffer.");
int *b_dest;
......@@ -857,9 +911,76 @@ void engine_redistribute(struct engine *e) {
}
swift_free("b_dest", b_dest);
/* Get destination of each sink-particle */
int *sink_counts;
if ((sink_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
error("Failed to allocate sink_counts temporary buffer.");
int *sink_dest;
if ((sink_dest = (int *)swift_malloc("sink_dest", sizeof(int) * nr_sinks)) ==
NULL)
error("Failed to allocate sink_dest temporary buffer.");
redist_data.counts = sink_counts;
redist_data.dest = sink_dest;
redist_data.base = (void *)sinks;
threadpool_map(&e->threadpool, engine_redistribute_dest_mapper_sink, sinks,
nr_sinks, sizeof(struct sink), threadpool_auto_chunk_size,
&redist_data);
/* Sort the particles according to their cell index. */
if (nr_sinks > 0)
space_sinks_sort(s->sinks, sink_dest, &sink_counts[nodeID * nr_nodes],
nr_nodes, 0);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the sink have been sorted correctly. */
for (size_t k = 0; k < nr_sinks; k++) {
const struct sink *sink = &s->sinks[k];
if (sink->time_bin == time_bin_inhibited)
error("Inhibited particle found after sorting!");
if (sink->time_bin == time_bin_not_created)
error("Inhibited particle found after sorting!");
/* New cell index */
const int new_cid =
cell_getid(s->cdim, sink->x[0] * s->iwidth[0],
sink->x[1] * s->iwidth[1], sink->x[2] * s->iwidth[2]);
/* New cell of this sink */
const struct cell *c = &s->cells_top[new_cid];
const int new_node = c->nodeID;
if (sink_dest[k] != new_node)
error("sink's new node index not matching sorted index.");
if (sink->x[0] < c->loc[0] || sink->x[0] > c->loc[0] + c->width[0] ||
sink->x[1] < c->loc[1] || sink->x[1] > c->loc[1] + c->width[1] ||
sink->x[2] < c->loc[2] || sink->x[2] > c->loc[2] + c->width[2])
error("sink not sorted into the right top-level cell!");
}
#endif
/* We need to re-link the gpart partners of sinks. */
if (nr_sinks > 0) {
struct savelink_mapper_data savelink_data;
savelink_data.nr_nodes = nr_nodes;
savelink_data.counts = sink_counts;
savelink_data.parts = (void *)sinks;
savelink_data.nodeID = nodeID;
threadpool_map(&e->threadpool, engine_redistribute_savelink_mapper_sink,
nodes, nr_nodes, sizeof(int), threadpool_auto_chunk_size,
&savelink_data);
}
swift_free("sink_dest", sink_dest);
/* Get destination of each g-particle */
int *g_counts;
if ((g_counts = (int *)calloc(sizeof(int), nr_nodes * nr_nodes)) == NULL)
if ((g_counts = (int *)calloc(nr_nodes * nr_nodes, sizeof(int))) == NULL)
error("Failed to allocate g_gcount temporary buffer.");
int *g_dest;
......@@ -932,22 +1053,30 @@ void engine_redistribute(struct engine *e) {
MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce bparticle transfer counts.");
/* Get all the sink_counts from all the nodes. */
if (MPI_Allreduce(MPI_IN_PLACE, sink_counts, nr_nodes * nr_nodes, MPI_INT,
MPI_SUM, MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to allreduce sink particle transfer counts.");
/* Report how many particles will be moved. */
if (e->verbose) {
if (e->nodeID == 0) {
size_t total = 0, g_total = 0, s_total = 0, b_total = 0;
size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0;
size_t total = 0, g_total = 0, s_total = 0, b_total = 0, sink_total = 0;
size_t unmoved = 0, g_unmoved = 0, s_unmoved = 0, b_unmoved = 0,
sink_unmoved = 0;
for (int p = 0, r = 0; p < nr_nodes; p++) {
for (int n = 0; n < nr_nodes; n++) {
total += counts[r];
g_total += g_counts[r];
s_total += s_counts[r];
b_total += b_counts[r];
sink_total += sink_counts[r];
if (p == n) {
unmoved += counts[r];
g_unmoved += g_counts[r];
s_unmoved += s_counts[r];
b_unmoved += b_counts[r];
sink_unmoved += sink_counts[r];
}
r++;
}
......@@ -967,14 +1096,19 @@ void engine_redistribute(struct engine *e) {
message("%ld of %ld (%.2f%%) of b-particles moved", b_total - b_unmoved,
b_total,
100.0 * (double)(b_total - b_unmoved) / (double)b_total);
if (sink_total > 0)
message(
"%ld of %ld (%.2f%%) of sink-particles moved",
sink_total - sink_unmoved, sink_total,
100.0 * (double)(sink_total - sink_unmoved) / (double)sink_total);
}
}
/* Now each node knows how many parts, sparts, bparts, and gparts will be
* transferred to every other node. Get the new numbers of particles for this
* node. */
/* Now each node knows how many parts, sparts, bparts, sinks and gparts will
* be transferred to every other node. Get the new numbers of particles for
* this node. */
size_t nr_parts_new = 0, nr_gparts_new = 0, nr_sparts_new = 0,
nr_bparts_new = 0;
nr_bparts_new = 0, nr_sinks_new = 0;
for (int k = 0; k < nr_nodes; k++)
nr_parts_new += counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
......@@ -983,6 +1117,8 @@ void engine_redistribute(struct engine *e) {
nr_sparts_new += s_counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
nr_bparts_new += b_counts[k * nr_nodes + nodeID];
for (int k = 0; k < nr_nodes; k++)
nr_sinks_new += sink_counts[k * nr_nodes + nodeID];
#ifdef WITH_CSDS
const int initial_redistribute = e->ti_current == 0;
......@@ -992,6 +1128,7 @@ void engine_redistribute(struct engine *e) {
size_t spart_offset = 0;
size_t gpart_offset = 0;
size_t bpart_offset = 0;
size_t sink_offset = 0;
for (int i = 0; i < nr_nodes; i++) {
const size_t c_ind = engine_rank * nr_nodes + i;
......@@ -1002,6 +1139,7 @@ void engine_redistribute(struct engine *e) {
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
sink_offset += sink_counts[c_ind];
continue;
}
......@@ -1023,11 +1161,17 @@ void engine_redistribute(struct engine *e) {
error("TODO");
}
/* Log the sinks */
if (sink_counts[c_ind] > 0) {
error("TODO");
}
/* Update the counters */
part_offset += counts[c_ind];
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
sink_offset += sink_counts[c_ind];
}
}
#endif
......@@ -1081,6 +1225,15 @@ void engine_redistribute(struct engine *e) {
s->nr_bparts = nr_bparts_new;
s->size_bparts = engine_redistribute_alloc_margin * nr_bparts_new;
/* Sink particles. */
new_parts = engine_do_redistribute(
"sinks", sink_counts, (char *)s->sinks, nr_sinks_new, sizeof(struct sink),
sink_align, sink_mpi_type, nr_nodes, nodeID, e->syncredist);
swift_free("sinks", s->sinks);
s->sinks = (struct sink *)new_parts;
s->nr_sinks = nr_sinks_new;
s->size_sinks = engine_redistribute_alloc_margin * nr_sinks_new;
/* All particles have now arrived. Time for some final operations on the
stuff we just received */
......@@ -1090,6 +1243,7 @@ void engine_redistribute(struct engine *e) {
size_t spart_offset = 0;
size_t gpart_offset = 0;
size_t bpart_offset = 0;
size_t sink_offset = 0;
for (int i = 0; i < nr_nodes; i++) {
const size_t c_ind = i * nr_nodes + engine_rank;
......@@ -1100,6 +1254,7 @@ void engine_redistribute(struct engine *e) {
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
sink_offset += sink_counts[c_ind];
continue;
}
......@@ -1121,11 +1276,17 @@ void engine_redistribute(struct engine *e) {
error("TODO");
}
/* Log the sinks */
if (sink_counts[c_ind] > 0) {
error("TODO");
}
/* Update the counters */
part_offset += counts[c_ind];
spart_offset += s_counts[c_ind];
gpart_offset += g_counts[c_ind];
bpart_offset += b_counts[c_ind];
sink_offset += sink_counts[c_ind];
}
}
#endif
......@@ -1139,6 +1300,7 @@ void engine_redistribute(struct engine *e) {
relink_data.g_counts = g_counts;
relink_data.s_counts = s_counts;
relink_data.b_counts = b_counts;
relink_data.sink_counts = sink_counts;
relink_data.nodeID = nodeID;
relink_data.nr_nodes = nr_nodes;
......@@ -1151,6 +1313,7 @@ void engine_redistribute(struct engine *e) {
free(g_counts);
free(s_counts);
free(b_counts);
free(sink_counts);
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that all parts are in the right place. */
......@@ -1186,6 +1349,15 @@ void engine_redistribute(struct engine *e) {
error("Received b-particle (%zu) that does not belong here (nodeID=%i).",
k, cells[cid].nodeID);
}
for (size_t k = 0; k < nr_sinks_new; k++) {
const int cid = cell_getid(s->cdim, s->sinks[k].x[0] * s->iwidth[0],
s->sinks[k].x[1] * s->iwidth[1],
s->sinks[k].x[2] * s->iwidth[2]);
if (cells[cid].nodeID != nodeID)
error(
"Received sink-particle (%zu) that does not belong here (nodeID=%i).",
k, cells[cid].nodeID);
}
/* Verify that the links are correct */
part_verify_links(s->parts, s->gparts, s->sinks, s->sparts, s->bparts,
......@@ -1200,10 +1372,11 @@ void engine_redistribute(struct engine *e) {
for (int k = 0; k < nr_cells; k++)
if (cells[k].nodeID == nodeID) my_cells += 1;
message(
"node %i now has %zu parts, %zu sparts, %zu bparts and %zu gparts in "
"node %i now has %zu parts, %zu sparts, %zu bparts, %zu sinks and %zu "
"gparts in "
"%i cells.",
nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_gparts_new,
my_cells);
nodeID, nr_parts_new, nr_sparts_new, nr_bparts_new, nr_sinks_new,
nr_gparts_new, my_cells);
}
/* Flag that we do not have any extra particles any more */
......@@ -1211,6 +1384,7 @@ void engine_redistribute(struct engine *e) {
s->nr_extra_gparts = 0;
s->nr_extra_sparts = 0;
s->nr_extra_bparts = 0;
s->nr_extra_sinks = 0;
/* Flag that a redistribute has taken place */
e->step_props |= engine_step_prop_redistribute;
......
......@@ -33,6 +33,7 @@
#include "random.h"
#include "rt.h"
#include "star_formation.h"
#include "tools.h"
#include "tracers.h"
const int particle_split_factor = 2;
......@@ -60,6 +61,7 @@ struct data_split {
long long offset_id;
long long *count_id;
swift_lock_type lock;
FILE *extra_split_logger;
};
/**
......@@ -172,19 +174,20 @@ void engine_split_gas_particle_split_mapper(void *restrict map_data, int count,
memcpy(&global_gparts[k_gparts], gp, sizeof(struct gpart));
}
/* Update splitting tree */
particle_splitting_update_binary_tree(&global_xparts[k_parts].split_data,
&xp->split_data);
/* Update the IDs. */
if (generate_random_ids) {
/* The gas IDs are always odd, so we multiply by two here to
* repsect the parity. */
* respect the parity. */
global_parts[k_parts].id += 2 * (long long)rand_r(&seedp);
} else {
global_parts[k_parts].id = offset_id + 2 * atomic_inc(count_id);
}
/* Update splitting tree */
particle_splitting_update_binary_tree(
&xp->split_data, &global_xparts[k_parts].split_data, p->id,
global_parts[k_parts].id, data->extra_split_logger, &data->lock);
/* Re-link everything */
if (with_gravity) {
global_parts[k_parts].gpart = &global_gparts[k_gparts];
......@@ -312,7 +315,10 @@ void engine_split_gas_particles(struct engine *e) {
/* Verify that nothing wrong happened with the IDs */
if (data_count.max_id > e->max_parts_id) {
error("Found a gas particle with an ID larger than the current max!");
error(
"Found a gas particle with an ID (%lld) larger than the current max "
"(%lld)!",
data_count.max_id, e->max_parts_id);
}
/* Be verbose about this. This is an important event */
......@@ -432,12 +438,22 @@ void engine_split_gas_particles(struct engine *e) {
size_t k_parts = s->nr_parts;
size_t k_gparts = s->nr_gparts;
FILE *extra_split_logger = NULL;
if (e->hydro_properties->log_extra_splits_in_file) {
char extra_split_logger_filename[256];
sprintf(extra_split_logger_filename, "splits/splits_%04d.txt", engine_rank);
extra_split_logger = fopen(extra_split_logger_filename, "a");
if (extra_split_logger == NULL) error("Error opening split logger file!");
}
/* Loop over the particles again to split them */
long long local_count_id = 0;
struct data_split data_split = {
e, mass_threshold, generate_random_ids, &k_parts,
&k_gparts, offset_id, &local_count_id, 0};
e, mass_threshold, generate_random_ids, &k_parts,
&k_gparts, offset_id, &local_count_id,
/*lock=*/0, extra_split_logger};
lock_init(&data_split.lock);
threadpool_map(&e->threadpool, engine_split_gas_particle_split_mapper,
s->parts, nr_parts_old, sizeof(struct part), 0, &data_split);
if (lock_destroy(&data_split.lock) != 0) error("Error destroying lock");
......@@ -459,7 +475,33 @@ void engine_split_gas_particles(struct engine *e) {
}
#endif
/* Close the logger file */
if (e->hydro_properties->log_extra_splits_in_file) fclose(extra_split_logger);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
void engine_init_split_gas_particles(struct engine *e) {
if (e->hydro_properties->log_extra_splits_in_file) {
/* Create the directory to host the logs */
if (engine_rank == 0) safe_checkdir("splits", /*create=*/1);
#ifdef WITH_MPI
MPI_Barrier(MPI_COMM_WORLD);
#endif
/* Create the logger files and add a header */
char extra_split_logger_filename[256];
sprintf(extra_split_logger_filename, "splits/splits_%04d.txt", engine_rank);
FILE *extra_split_logger = fopen(extra_split_logger_filename, "w");
fprintf(extra_split_logger, "# %12s %20s %20s %20s %20s\n", "Step", "ID",
"Progenitor", "Count", "Tree");
/* Close everything for now */
fclose(extra_split_logger);
}
}
......@@ -33,6 +33,13 @@
/* Local headers. */
#include "proxy.h"
#ifdef WITH_MPI
/* Number of particle types to wait for after launching the proxies. We have
parts, xparts, gparts, sparts, bparts and sinks to exchange, hence 6 types.
*/
#define MPI_REQUEST_NUMBER_PARTICLE_TYPES 6
#endif
/**
* @brief Exchange straying particles with other nodes.
*
......@@ -57,18 +64,22 @@
* @param ind_bpart The foreign #cell ID of each bpart.
* @param Nbpart The number of stray bparts, contains the number of bparts
* received on return.
* @param offset_sinks The index in the sinks array as of which the foreign
* parts reside (i.e. the current number of local #sink).
* @param ind_sink The foreign #cell ID of each sink.
* @param Nsink The number of stray sinks, contains the number of sinks
* received on return.
*
* Note that this function does not mess-up the linkage between parts and
* gparts, i.e. the received particles have correct linkeage.
*/
void engine_exchange_strays(struct engine *e, const size_t offset_parts,
const int *restrict ind_part, size_t *Npart,
const size_t offset_gparts,
const int *restrict ind_gpart, size_t *Ngpart,
const size_t offset_sparts,
const int *restrict ind_spart, size_t *Nspart,
const size_t offset_bparts,
const int *restrict ind_bpart, size_t *Nbpart) {
void engine_exchange_strays(
struct engine *e, const size_t offset_parts, const int *restrict ind_part,
size_t *Npart, const size_t offset_gparts, const int *restrict ind_gpart,
size_t *Ngpart, const size_t offset_sparts, const int *restrict ind_spart,
size_t *Nspart, const size_t offset_bparts, const int *restrict ind_bpart,
size_t *Nbpart, const size_t offset_sinks, const int *restrict ind_sink,
size_t *Nsink) {
#ifdef WITH_MPI
struct space *s = e->s;
......@@ -80,6 +91,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
e->proxies[k].nr_gparts_out = 0;
e->proxies[k].nr_sparts_out = 0;
e->proxies[k].nr_bparts_out = 0;
e->proxies[k].nr_sinks_out = 0;
}
/* Put the parts into the corresponding proxies. */
......@@ -204,6 +216,46 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
/* Load the bpart into the proxy */
proxy_bparts_load(&e->proxies[pid], &s->bparts[offset_bparts + k], 1);
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
error("Not yet implemented.");
}
#endif
}
/* Put the sinks into the corresponding proxies. */
for (size_t k = 0; k < *Nsink; k++) {
/* Ignore the particles we want to get rid of (inhibited, ...). */
if (ind_sink[k] == -1) continue;
/* Get the target node and proxy ID. */
const int node_id = e->s->cells_top[ind_sink[k]].nodeID;
if (node_id < 0 || node_id >= e->nr_nodes)
error("Bad node ID %i.", node_id);
const int pid = e->proxy_ind[node_id];
if (pid < 0) {
error(
"Do not have a proxy for the requested nodeID %i for part with "
"id=%lld, x=[%e,%e,%e].",
node_id, s->sinks[offset_sinks + k].id,
s->sinks[offset_sinks + k].x[0], s->sinks[offset_sinks + k].x[1],
s->sinks[offset_sinks + k].x[2]);
}
/* Re-link the associated gpart with the buffer offset of the sink. */
if (s->sinks[offset_sinks + k].gpart != NULL) {
s->sinks[offset_sinks + k].gpart->id_or_neg_offset =
-e->proxies[pid].nr_sinks_out;
}
#ifdef SWIFT_DEBUG_CHECKS
if (s->sinks[offset_sinks + k].time_bin == time_bin_inhibited)
error("Attempting to exchange an inhibited particle");
#endif
/* Load the sink into the proxy */
proxy_sinks_load(&e->proxies[pid], &s->sinks[offset_sinks + k], 1);
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
error("Not yet implemented.");
......@@ -252,8 +304,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
}
/* Launch the proxies. */
MPI_Request reqs_in[5 * engine_maxproxies];
MPI_Request reqs_out[5 * engine_maxproxies];
MPI_Request reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
MPI_Request reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * engine_maxproxies];
for (int k = 0; k < e->nr_proxies; k++) {
proxy_parts_exchange_first(&e->proxies[k]);
reqs_in[k] = e->proxies[k].req_parts_count_in;
......@@ -281,18 +333,21 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
int count_gparts_in = 0;
int count_sparts_in = 0;
int count_bparts_in = 0;
int count_sinks_in = 0;
for (int k = 0; k < e->nr_proxies; k++) {
count_parts_in += e->proxies[k].nr_parts_in;
count_gparts_in += e->proxies[k].nr_gparts_in;
count_sparts_in += e->proxies[k].nr_sparts_in;
count_bparts_in += e->proxies[k].nr_bparts_in;
count_sinks_in += e->proxies[k].nr_sinks_in;
}
if (e->verbose) {
message(
"sent out %zu/%zu/%zu/%zu parts/gparts/sparts/bparts, got %i/%i/%i/%i "
"sent out %zu/%zu/%zu/%zu/%zu parts/gparts/sparts/bparts/sinks, got "
"%i/%i/%i/%i/%i "
"back.",
*Npart, *Ngpart, *Nspart, *Nbpart, count_parts_in, count_gparts_in,
count_sparts_in, count_bparts_in);
*Npart, *Ngpart, *Nspart, *Nbpart, *Nsink, count_parts_in,
count_gparts_in, count_sparts_in, count_bparts_in, count_sinks_in);
}
/* Reallocate the particle arrays if necessary */
......@@ -356,6 +411,24 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
}
}
if (offset_sinks + count_sinks_in > s->size_sinks) {
s->size_sinks = (offset_sinks + count_sinks_in) * engine_parts_size_grow;
struct sink *sinks_new = NULL;
if (swift_memalign("sinks", (void **)&sinks_new, sink_align,
sizeof(struct sink) * s->size_sinks) != 0)
error("Failed to allocate new sink data.");
memcpy(sinks_new, s->sinks, sizeof(struct sink) * offset_sinks);
swift_free("sinks", s->sinks);
s->sinks = sinks_new;
/* Reset the links */
for (size_t k = 0; k < offset_sinks; k++) {
if (s->sinks[k].gpart != NULL) {
s->sinks[k].gpart->id_or_neg_offset = -k;
}
}
}
if (offset_gparts + count_gparts_in > s->size_gparts) {
s->size_gparts = (offset_gparts + count_gparts_in) * engine_parts_size_grow;
struct gpart *gparts_new = NULL;
......@@ -374,6 +447,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_black_hole) {
s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_sink) {
s->sinks[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
}
}
......@@ -382,82 +457,113 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
int nr_in = 0, nr_out = 0;
for (int k = 0; k < e->nr_proxies; k++) {
if (e->proxies[k].nr_parts_in > 0) {
reqs_in[5 * k] = e->proxies[k].req_parts_in;
reqs_in[5 * k + 1] = e->proxies[k].req_xparts_in;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
e->proxies[k].req_parts_in;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
e->proxies[k].req_xparts_in;
nr_in += 2;
} else {
reqs_in[5 * k] = reqs_in[5 * k + 1] = MPI_REQUEST_NULL;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_gparts_in > 0) {
reqs_in[5 * k + 2] = e->proxies[k].req_gparts_in;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
e->proxies[k].req_gparts_in;
nr_in += 1;
} else {
reqs_in[5 * k + 2] = MPI_REQUEST_NULL;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sparts_in > 0) {
reqs_in[5 * k + 3] = e->proxies[k].req_sparts_in;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
e->proxies[k].req_sparts_in;
nr_in += 1;
} else {
reqs_in[5 * k + 3] = MPI_REQUEST_NULL;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_bparts_in > 0) {
reqs_in[5 * k + 4] = e->proxies[k].req_bparts_in;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
e->proxies[k].req_bparts_in;
nr_in += 1;
} else {
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sinks_in > 0) {
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
e->proxies[k].req_sinks_in;
nr_in += 1;
} else {
reqs_in[5 * k + 4] = MPI_REQUEST_NULL;
reqs_in[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_parts_out > 0) {
reqs_out[5 * k] = e->proxies[k].req_parts_out;
reqs_out[5 * k + 1] = e->proxies[k].req_xparts_out;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
e->proxies[k].req_parts_out;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
e->proxies[k].req_xparts_out;
nr_out += 2;
} else {
reqs_out[5 * k] = reqs_out[5 * k + 1] = MPI_REQUEST_NULL;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k] =
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 1] =
MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_gparts_out > 0) {
reqs_out[5 * k + 2] = e->proxies[k].req_gparts_out;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] =
e->proxies[k].req_gparts_out;
nr_out += 1;
} else {
reqs_out[5 * k + 2] = MPI_REQUEST_NULL;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 2] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sparts_out > 0) {
reqs_out[5 * k + 3] = e->proxies[k].req_sparts_out;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] =
e->proxies[k].req_sparts_out;
nr_out += 1;
} else {
reqs_out[5 * k + 3] = MPI_REQUEST_NULL;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 3] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_bparts_out > 0) {
reqs_out[5 * k + 4] = e->proxies[k].req_bparts_out;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] =
e->proxies[k].req_bparts_out;
nr_out += 1;
} else {
reqs_out[5 * k + 4] = MPI_REQUEST_NULL;
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 4] = MPI_REQUEST_NULL;
}
if (e->proxies[k].nr_sinks_out > 0) {
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] =
e->proxies[k].req_sinks_out;
nr_out += 1;
} else {
reqs_out[MPI_REQUEST_NUMBER_PARTICLE_TYPES * k + 5] = MPI_REQUEST_NULL;
}
}
/* Wait for each part array to come in and collect the new
parts from the proxies. */
int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0;
int count_parts = 0, count_gparts = 0, count_sparts = 0, count_bparts = 0,
count_sinks = 0;
for (int k = 0; k < nr_in; k++) {
int err, pid;
if ((err = MPI_Waitany(5 * e->nr_proxies, reqs_in, &pid,
MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
if ((err = MPI_Waitany(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies,
reqs_in, &pid, MPI_STATUS_IGNORE)) != MPI_SUCCESS) {
char buff[MPI_MAX_ERROR_STRING];
int res;
MPI_Error_string(err, buff, &res);
error("MPI_Waitany failed (%s).", buff);
}
if (pid == MPI_UNDEFINED) break;
// message( "request from proxy %i has arrived." , pid / 5 );
pid = 5 * (pid / 5);
// message( "request from proxy %i has arrived." , pid /
// MPI_REQUEST_NUMBER_PARTICLE_TYPES );
pid = MPI_REQUEST_NUMBER_PARTICLE_TYPES *
(pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES);
/* If all the requests for a given proxy have arrived... */
if (reqs_in[pid + 0] == MPI_REQUEST_NULL &&
reqs_in[pid + 1] == MPI_REQUEST_NULL &&
reqs_in[pid + 2] == MPI_REQUEST_NULL &&
reqs_in[pid + 3] == MPI_REQUEST_NULL &&
reqs_in[pid + 4] == MPI_REQUEST_NULL) {
reqs_in[pid + 4] == MPI_REQUEST_NULL &&
reqs_in[pid + 5] == MPI_REQUEST_NULL) {
/* Copy the particle data to the part/xpart/gpart arrays. */
struct proxy *prox = &e->proxies[pid / 5];
struct proxy *prox = &e->proxies[pid / MPI_REQUEST_NUMBER_PARTICLE_TYPES];
memcpy(&s->parts[offset_parts + count_parts], prox->parts_in,
sizeof(struct part) * prox->nr_parts_in);
memcpy(&s->xparts[offset_parts + count_parts], prox->xparts_in,
......@@ -468,6 +574,8 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
sizeof(struct spart) * prox->nr_sparts_in);
memcpy(&s->bparts[offset_bparts + count_bparts], prox->bparts_in,
sizeof(struct bpart) * prox->nr_bparts_in);
memcpy(&s->sinks[offset_sinks + count_sinks], prox->sinks_in,
sizeof(struct sink) * prox->nr_sinks_in);
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
......@@ -495,6 +603,12 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
if (prox->nr_bparts_in > 0) {
error("TODO");
}
/* Log the sinks */
if (prox->nr_sinks_in > 0) {
/* Not implemented yet */
error("TODO");
}
}
#endif
/* for (int k = offset; k < offset + count; k++)
......@@ -522,6 +636,11 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
&s->bparts[offset_bparts + count_bparts - gp->id_or_neg_offset];
gp->id_or_neg_offset = s->bparts - bp;
bp->gpart = gp;
} else if (gp->type == swift_type_sink) {
struct sink *sink =
&s->sinks[offset_sinks + count_sinks - gp->id_or_neg_offset];
gp->id_or_neg_offset = s->sinks - sink;
sink->gpart = gp;
}
}
......@@ -530,13 +649,14 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
count_gparts += prox->nr_gparts_in;
count_sparts += prox->nr_sparts_in;
count_bparts += prox->nr_bparts_in;
count_sinks += prox->nr_sinks_in;
}
}
/* Wait for all the sends to have finished too. */
if (nr_out > 0)
if (MPI_Waitall(5 * e->nr_proxies, reqs_out, MPI_STATUSES_IGNORE) !=
MPI_SUCCESS)
if (MPI_Waitall(MPI_REQUEST_NUMBER_PARTICLE_TYPES * e->nr_proxies, reqs_out,
MPI_STATUSES_IGNORE) != MPI_SUCCESS)
error("MPI_Waitall on sends failed.");
/* Free the proxy memory */
......@@ -553,6 +673,7 @@ void engine_exchange_strays(struct engine *e, const size_t offset_parts,
*Ngpart = count_gparts;
*Nspart = count_sparts;
*Nbpart = count_bparts;
*Nsink = count_sinks;
#else
error("SWIFT was not compiled with MPI support.");
......
......@@ -160,12 +160,6 @@ static void engine_do_unskip_black_holes(struct cell *c, struct engine *e) {
/* Early abort (are we below the level where tasks are)? */
if (!cell_get_flag(c, cell_flag_has_tasks)) return;
/* Ignore empty cells. */
if (c->black_holes.count == 0) return;
/* Skip inactive cells. */
if (!cell_is_active_black_holes(c, e)) return;
/* Recurse */
if (c->split) {
for (int k = 0; k < 8; k++) {
......@@ -428,12 +422,6 @@ void engine_unskip(struct engine *e) {
memswap(&local_cells[k], &local_cells[num_active_cells], sizeof(int));
num_active_cells += 1;
}
/* Activate the top-level timestep exchange */
#ifdef WITH_MPI
scheduler_activate_all_subtype(&e->sched, c->mpi.send, task_subtype_tend);
scheduler_activate_all_subtype(&e->sched, c->mpi.recv, task_subtype_tend);
#endif
}
/* What kind of tasks do we have? */
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#include <config.h>
#ifdef HAVE_HDF5
#include <hdf5.h>
#endif
/* Local includes */
#include "cosmology.h"
#include "entropy_floor.h"
#include "hydro.h"
#include "parser.h"
/**
* @brief Compute the pressure from the entropy floor at a given density
*
* @param rho_phys The physical density (internal units).
* @param rho_com The comoving density (internal units).
* @param cosmo The cosmological model.
* @param props The properties of the entropy floor.
*/
float entropy_floor_gas_pressure(const float rho_phys, const float rho_com,
const struct cosmology *cosmo,
const struct entropy_floor_properties *props) {
/* Mean baryon density in co-moving internal units for over-density condition
* (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
* making the over-density condition a no-op) */
const float rho_crit_0 = cosmo->critical_density_0;
const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
/* Physical pressure */
float pressure = 0.f;
/* Are we in the regime of the Jeans equation of state? */
if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
(rho_phys >= props->Jeans_density_threshold)) {
const float pressure_Jeans =
props->Jeans_pressure_norm *
powf(rho_phys * props->Jeans_density_threshold_inv,
props->Jeans_gamma_effective);
pressure = max(pressure, pressure_Jeans);
}
/* Are we in the regime of the Cool equation of state? */
if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
(rho_phys >= props->Cool_density_threshold)) {
const float pressure_Cool =
props->Cool_pressure_norm *
powf(rho_phys * props->Cool_density_threshold_inv,
props->Cool_gamma_effective);
pressure = max(pressure, pressure_Cool);
}
return pressure;
}
/**
* @brief Compute the entropy floor of a given #part.
*
* Note that the particle is not updated!!
*
* @param p The #part.
* @param cosmo The cosmological model.
* @param props The properties of the entropy floor.
*/
float entropy_floor(const struct part *p, const struct cosmology *cosmo,
const struct entropy_floor_properties *props) {
/* Comoving density in internal units */
const float rho_com = hydro_get_comoving_density(p);
/* Physical density in internal units */
const float rho_phys = hydro_get_physical_density(p, cosmo);
const float pressure =
entropy_floor_gas_pressure(rho_phys, rho_com, cosmo, props);
/* Convert to an entropy.
* (Recall that the entropy is the same in co-moving and physical frames) */
return gas_entropy_from_pressure(rho_phys, pressure);
}
/**
* @brief Compute the temperature from the entropy floor at a given density
*
* This is the temperature exactly corresponding to the imposed EoS shape.
* It only matches the entropy returned by the entropy_floor() function
* for a neutral gas with primoridal abundance.
*
* @param rho_phys The physical density (internal units).
* @param rho_com The comoving density (internal units).
* @param cosmo The cosmological model.
* @param props The properties of the entropy floor.
*/
float entropy_floor_gas_temperature(
const float rho_phys, const float rho_com, const struct cosmology *cosmo,
const struct entropy_floor_properties *props) {
/* Mean baryon density in co-moving internal units for over-density condition
* (Recall cosmo->critical_density_0 is 0 in a non-cosmological run,
* making the over-density condition a no-op) */
const float rho_crit_0 = cosmo->critical_density_0;
const float rho_crit_baryon = cosmo->Omega_b * rho_crit_0;
/* Physical */
float temperature = 0.f;
/* Are we in the regime of the Jeans equation of state? */
if ((rho_com >= rho_crit_baryon * props->Jeans_over_density_threshold) &&
(rho_phys >= props->Jeans_density_threshold)) {
const float jeans_slope = props->Jeans_gamma_effective - 1.f;
const float temperature_Jeans =
props->Jeans_temperature_norm *
pow(rho_phys * props->Jeans_density_threshold_inv, jeans_slope);
temperature = max(temperature, temperature_Jeans);
}
/* Are we in the regime of the Cool equation of state? */
if ((rho_com >= rho_crit_baryon * props->Cool_over_density_threshold) &&
(rho_phys >= props->Cool_density_threshold)) {
const float cool_slope = props->Cool_gamma_effective - 1.f;
const float temperature_Cool =
props->Cool_temperature_norm *
pow(rho_phys * props->Cool_density_threshold_inv, cool_slope);
temperature = max(temperature, temperature_Cool);
}
return temperature;
}
/**
* @brief Compute the temperature from the entropy floor for a given #part
*
* Calculate the EoS temperature, the particle is not updated.
* This is the temperature exactly corresponding to the imposed EoS shape.
* It only matches the entropy returned by the entropy_floor() function
* for a neutral gas with primoridal abundance.
*
* @param p The #part.
* @param cosmo The cosmological model.
* @param props The properties of the entropy floor.
*/
float entropy_floor_temperature(const struct part *p,
const struct cosmology *cosmo,
const struct entropy_floor_properties *props) {
/* Comoving density in internal units */
const float rho_com = hydro_get_comoving_density(p);
/* Physical density in internal units */
const float rho_phys = hydro_get_physical_density(p, cosmo);
return entropy_floor_gas_temperature(rho_phys, rho_com, cosmo, props);
}
/**
* @brief Initialise the entropy floor by reading the parameters and converting
* to internal units.
*
* The input temperatures and number densities are converted to entropy and
* density assuming a neutral gas of primoridal abundance.
*
* @param params The YAML parameter file.
* @param us The system of units used internally.
* @param phys_const The physical constants.
* @param hydro_props The propoerties of the hydro scheme.
* @param props The entropy floor properties to fill.
*/
void entropy_floor_init(struct entropy_floor_properties *props,
const struct phys_const *phys_const,
const struct unit_system *us,
const struct hydro_props *hydro_props,
struct swift_params *params) {
/* Read the parameters in the units they are set */
props->Jeans_density_threshold_H_p_cm3 = parser_get_param_float(
params, "EAGLEEntropyFloor:Jeans_density_threshold_H_p_cm3");
props->Jeans_over_density_threshold = parser_get_param_float(
params, "EAGLEEntropyFloor:Jeans_over_density_threshold");
props->Jeans_temperature_norm_K = parser_get_param_float(
params, "EAGLEEntropyFloor:Jeans_temperature_norm_K");
props->Jeans_gamma_effective =
parser_get_param_float(params, "EAGLEEntropyFloor:Jeans_gamma_effective");
props->Cool_density_threshold_H_p_cm3 = parser_get_param_float(
params, "EAGLEEntropyFloor:Cool_density_threshold_H_p_cm3");
props->Cool_over_density_threshold = parser_get_param_float(
params, "EAGLEEntropyFloor:Cool_over_density_threshold");
props->Cool_temperature_norm_K = parser_get_param_float(
params, "EAGLEEntropyFloor:Cool_temperature_norm_K");
props->Cool_gamma_effective =
parser_get_param_float(params, "EAGLEEntropyFloor:Cool_gamma_effective");
/* Cross-check that the input makes sense */
if (props->Cool_density_threshold_H_p_cm3 >=
props->Jeans_density_threshold_H_p_cm3) {
error(
"Invalid values for the entrop floor density thresholds. The 'Jeans' "
"threshold (%e cm^-3) should be at a higher density than the 'Cool' "
"threshold (%e cm^-3)",
props->Jeans_density_threshold_H_p_cm3,
props->Cool_density_threshold_H_p_cm3);
}
/* Initial Hydrogen abundance (mass fraction) */
const double X_H = hydro_props->hydrogen_mass_fraction;
/* Now convert to internal units assuming primodial Hydrogen abundance */
props->Jeans_temperature_norm =
props->Jeans_temperature_norm_K /
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
props->Jeans_density_threshold =
props->Jeans_density_threshold_H_p_cm3 /
units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
phys_const->const_proton_mass / X_H;
props->Cool_temperature_norm =
props->Cool_temperature_norm_K /
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
props->Cool_density_threshold =
props->Cool_density_threshold_H_p_cm3 /
units_cgs_conversion_factor(us, UNIT_CONV_NUMBER_DENSITY) *
phys_const->const_proton_mass / X_H;
/* We assume neutral gas */
const float mean_molecular_weight = hydro_props->mu_neutral;
/* Get the common terms */
props->Jeans_density_threshold_inv = 1.f / props->Jeans_density_threshold;
props->Cool_density_threshold_inv = 1.f / props->Cool_density_threshold;
/* P_norm = (k_B * T) / (m_p * mu) * rho_threshold */
props->Jeans_pressure_norm =
((phys_const->const_boltzmann_k * props->Jeans_temperature_norm) /
(phys_const->const_proton_mass * mean_molecular_weight)) *
props->Jeans_density_threshold;
props->Cool_pressure_norm =
((phys_const->const_boltzmann_k * props->Cool_temperature_norm) /
(phys_const->const_proton_mass * mean_molecular_weight)) *
props->Cool_density_threshold;
}
/**
* @brief Print the properties of the entropy floor to stdout.
*
* @param props The entropy floor properties.
*/
void entropy_floor_print(const struct entropy_floor_properties *props) {
message("Entropy floor is 'EAGLE' with:");
message("Jeans limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
props->Jeans_gamma_effective, props->Jeans_density_threshold,
props->Jeans_density_threshold_H_p_cm3,
props->Jeans_temperature_norm);
message(" Cool limiter with slope n=%.3f at rho=%e (%e H/cm^3) and T=%.1f K",
props->Cool_gamma_effective, props->Cool_density_threshold,
props->Cool_density_threshold_H_p_cm3, props->Cool_temperature_norm);
}
#ifdef HAVE_HDF5
/**
* @brief Writes the current model of entropy floor to the file
* @param h_grp The HDF5 group in which to write
*/
void entropy_floor_write_flavour(hid_t h_grp) {
io_write_attribute_s(h_grp, "Entropy floor", "EAGLE");
}
#endif
/**
* @brief Write an entropy floor struct to the given FILE as a stream of bytes.
*
* @param props the struct
* @param stream the file stream
*/
void entropy_floor_struct_dump(const struct entropy_floor_properties *props,
FILE *stream) {
restart_write_blocks((void *)props, sizeof(struct entropy_floor_properties),
1, stream, "entropy floor", "entropy floor properties");
}
/**
* @brief Restore a entropy floor struct from the given FILE as a stream of
* bytes.
*
* @param props the struct
* @param stream the file stream
*/
void entropy_floor_struct_restore(struct entropy_floor_properties *props,
FILE *stream) {
restart_read_blocks((void *)props, sizeof(struct entropy_floor_properties), 1,
stream, NULL, "entropy floor properties");
}