Commit cfd9a324 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'feedback' into 'master'

Star density

See merge request !581
parents d9fd1872 327eecab
......@@ -63,12 +63,18 @@ tests/brute_force_perturbed.dat
tests/swift_dopair_perturbed.dat
tests/test27cells
tests/test27cells_subset
tests/test27cellsStars
tests/test27cellsStars_subset
tests/testPeriodicBC
tests/test125cells
tests/brute_force_27_standard.dat
tests/swift_dopair_27_standard.dat
tests/brute_force_27_perturbed.dat
tests/swift_dopair_27_perturbed.dat
tests/star_brute_force_27_standard.dat
tests/swift_star_dopair_27_standard.dat
tests/star_brute_force_27_perturbed.dat
tests/swift_star_dopair_27_perturbed.dat
tests/brute_force_125_standard.dat
tests/swift_dopair_125_standard.dat
tests/brute_force_125_perturbed.dat
......@@ -106,6 +112,8 @@ tests/testPeriodicBC.sh
tests/testPeriodicBCPerturbed.sh
tests/test27cells.sh
tests/test27cellsPerturbed.sh
tests/test27cellsStars.sh
tests/test27cellsStarsPerturbed.sh
tests/test125cells.sh
tests/test125cellsPerturbed.sh
tests/testParser.sh
......
......@@ -1045,11 +1045,15 @@ case "$with_subgrid" in
with_subgrid_cooling=grackle
with_subgrid_chemistry=GEAR
with_subgrid_hydro=gadget2
with_subgrid_stars=GEAR
with_subgrid_feedback=thermal
;;
EAGLE)
with_subgrid_cooling=EAGLE
with_subgrid_chemistry=EAGLE
with_subgrid_hydro=gadget2
with_subgrid_stars=none
with_subgrid_feedback=none
;;
*)
AC_MSG_ERROR([Unknown subgrid choice: $with_subgrid])
......@@ -1129,6 +1133,25 @@ case "$with_hydro" in
;;
esac
# Check if debugging interactions stars is switched on.
AC_ARG_ENABLE([debug-interactions-stars],
[AS_HELP_STRING([--enable-debug-interactions-stars],
[Activate interaction debugging for stars, logging a maximum of @<:@N@:>@ neighbours. Defaults to 256 if no value set.]
)],
[enable_debug_interactions_stars="$enableval"],
[enable_debug_interactions_stars="no"]
)
if test "$enable_debug_interactions_stars" != "no"; then
AC_DEFINE([DEBUG_INTERACTIONS_STARS],1,[Enable interaction debugging for stars])
if test "$enable_debug_interactions_stars" == "yes"; then
AC_DEFINE([MAX_NUM_OF_NEIGHBOURS_STARS],256,[The maximum number of particle neighbours to be logged for stars])
[enable_debug_interactions_stars="yes (Logging up to 256 neighbours)"]
else
AC_DEFINE_UNQUOTED([MAX_NUM_OF_NEIGHBOURS_STARS], [$enableval] ,[The maximum number of particle neighbours to be logged for stars])
[enable_debug_interactions_stars="yes (Logging up to $enableval neighbours)"]
fi
fi
# Check if debugging interactions is switched on.
AC_ARG_ENABLE([debug-interactions],
[AS_HELP_STRING([--enable-debug-interactions],
......@@ -1152,6 +1175,7 @@ if test "$enable_debug_interactions" != "no"; then
fi
fi
# SPH Kernel function
AC_ARG_WITH([kernel],
[AS_HELP_STRING([--with-kernel=<kernel>],
......@@ -1401,6 +1425,64 @@ case "$with_chemistry" in
;;
esac
# Stellar model.
AC_ARG_WITH([stars],
[AS_HELP_STRING([--with-stars=<model>],
[Stellar model to use @<:@none, GEAR, debug default: none@:>@]
)],
[with_stars="$withval"],
[with_stars="none"]
)
if test "$with_subgrid" != "none"; then
if test "$with_stars" != "none"; then
AC_MSG_ERROR([Cannot provide with-subgrid and with-stars together])
else
with_stars="$with_subgrid_stars"
fi
fi
case "$with_stars" in
GEAR)
AC_DEFINE([STARS_GEAR], [1], [GEAR stellar model])
;;
none)
;;
*)
AC_MSG_ERROR([Unknown stellar model: $with_stars])
;;
esac
# Feedback model
AC_ARG_WITH([feedback],
[AS_HELP_STRING([--with-feedback=<model>],
[Feedback model to use @<:@none, thermal, debug default: none@:>@]
)],
[with_feedback="$withval"],
[with_feedback="none"]
)
if test "$with_subgrid" != "none"; then
if test "$with_feedback" != "none"; then
AC_MSG_ERROR([Cannot provide with-subgrid and with-feedback together])
else
with_feedback="$with_subgrid_feedback"
fi
fi
case "$with_feedback" in
thermal)
AC_DEFINE([FEEDBACK_THERMAL], [1], [Thermal Blastwave])
;;
none)
;;
*)
AC_MSG_ERROR([Unknown feedback model: $with_feedback])
;;
esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
......@@ -1460,6 +1542,8 @@ AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh])
AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh])
AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
AC_CONFIG_FILES([tests/test27cellsPerturbed.sh], [chmod +x tests/test27cellsPerturbed.sh])
AC_CONFIG_FILES([tests/test27cellsStars.sh], [chmod +x tests/test27cellsStars.sh])
AC_CONFIG_FILES([tests/test27cellsStarsPerturbed.sh], [chmod +x tests/test27cellsStarsPerturbed.sh])
AC_CONFIG_FILES([tests/test125cells.sh], [chmod +x tests/test125cells.sh])
AC_CONFIG_FILES([tests/test125cellsPerturbed.sh], [chmod +x tests/test125cellsPerturbed.sh])
AC_CONFIG_FILES([tests/testPeriodicBC.sh], [chmod +x tests/testPeriodicBC.sh])
......@@ -1521,6 +1605,8 @@ AC_MSG_RESULT([
Cooling function : $with_cooling
Chemistry : $with_chemistry
Stellar model : $with_stars
Feedback model : $with_feedback
Individual timers : $enable_timers
Task debugging : $enable_task_debugging
......
.. Task
Loic Hausammann 17th July 2018
.. _task_adding_your_own:
.. highlight:: c
Adding a Task
=============
First you will need to understand the dependencies between tasks
using the file ``dependency_graph.dot`` generated by swift at the beginning of any simulation and then decide where it will fit (see :ref:`task`).
For the next paragraphs, let's assume that we want to implement the existing task ``cooling``.
Adding it to the Task List
--------------------------
First you will need to add it to the task list situated in ``task.h`` and ``task.c``.
In ``task.h``, you need to provide an additional entry to the enum ``task_types`` (e.g. ``task_type_cooling``).
The last entry ``task_type_count`` should always stay at the end as it is a counter of the number of elements.
For example::
enum task_types {
task_type_none = 0,
task_type_sort,
task_type_self,
task_type_pair,
task_type_sub_self,
task_type_sub_pair,
task_type_ghost_in,
task_type_ghost,
task_type_ghost_out,
task_type_extra_ghost,
task_type_drift_part,
task_type_end_force,
task_type_kick1,
task_type_kick2,
task_type_timestep,
task_type_send,
task_type_recv,
task_type_cooling,
task_type_count
} __attribute__((packed));
In ``task.c``, you will find an array containing the name of each task and need to add your own (e.g. ``cooling``).
Be careful with the order that should be the same than in the previous list.
For example::
/* Task type names. */
const char *taskID_names[task_type_count] = {
"none", "sort", "self", "pair", "sub_self",
"sub_pair", "ghost_in", "ghost", "ghost_out",
"extra_ghost", "drift_part", "end_force", "kick1",
"kick2", "timestep", "send", "recv",
"cooling"};
Adding it to the Cells
----------------------
Each cell contains a list to its tasks and therefore you need to provide a link for it.
In ``cell.h``, add a pointer to a task in the structure.
In order to stay clean, please put the new task in the same group than the other tasks.
For example::
struct cell {
/* Lot of stuff before. */
/*! Task for the cooling */
struct task *cooling;
/*! The second kick task */
struct task *kick2;
/* Lot of stuff after */
}
Adding a new Timer
------------------
As SWIFT is HPC oriented, any new task need to be optimized.
It cannot be done without timing the function.
In ``timers.h``, you will find an enum that contains all the tasks.
You will need to add yours inside it.
For example::
enum {
timer_none = 0,
timer_prepare,
timer_init,
timer_drift_part,
timer_drift_gpart,
timer_kick1,
timer_kick2,
timer_timestep,
timer_endforce,
timer_dosort,
timer_doself_density,
timer_doself_gradient,
timer_doself_force,
timer_dopair_density,
timer_dopair_gradient,
timer_dopair_force,
timer_dosub_self_density,
timer_dosub_self_gradient,
timer_dosub_self_force,
timer_dosub_pair_density,
timer_dosub_pair_gradient,
timer_dosub_pair_force,
timer_doself_subset,
timer_dopair_subset,
timer_dopair_subset_naive,
timer_dosub_subset,
timer_do_ghost,
timer_do_extra_ghost,
timer_dorecv_part,
timer_do_cooling,
timer_gettask,
timer_qget,
timer_qsteal,
timer_locktree,
timer_runners,
timer_step,
timer_cooling,
timer_count,
};
As for ``task.h``,
you will need to give a name to your timer in ``timers.c``::
const char* timers_names[timer_count] = {
"none",
"prepare",
"init",
"drift_part",
"kick1",
"kick2",
"timestep",
"endforce",
"dosort",
"doself_density",
"doself_gradient",
"doself_force",
"dopair_density",
"dopair_gradient",
"dopair_force",
"dosub_self_density",
"dosub_self_gradient",
"dosub_self_force",
"dosub_pair_density",
"dosub_pair_gradient",
"dosub_pair_force",
"doself_subset",
"dopair_subset",
"dopair_subset_naive",
"dosub_subset",
"do_ghost",
"do_extra_ghost",
"dorecv_part",
"gettask",
"qget",
"qsteal",
"locktree",
"runners",
"step",
"cooling",
};
You can now easily time
your functions by using::
TIMER_TIC;
/* Your complicated functions */
if (timer) TIMER_TOC(timer_cooling);
Adding your Task to the System
------------------------------
Now the tricky part happens.
SWIFT is able to deal automatically with the conflicts between tasks, but unfortunately cannot understand the dependencies.
To implement your new task in the task system, you will need to modify a few functions in ``engine.c``.
First, you will need to add mainly two functions: ``scheduler_addtask`` and ``scheduler_addunlocks`` in the ``engine_make_hierarchical_tasks_*`` functions (depending on the type of task you implement, you will need to write it to a different function).
In ``engine_make_hierarchical_tasks_hydro``,
we add the task through the following call::
/* Add the cooling task. */
c->cooling =
scheduler_addtask(s, task_type_cooling, task_subtype_none, 0,
0, c, NULL);
As the ``cooling`` cannot be done before the end of the force computation
and the second kick cannot be done before the cooling::
scheduler_addunlock(s, c->super->end_force, c->cooling);
scheduler_addunlock(s, c->cooling, c->super->kick2);
The next step is to activate your task
in ``engine_marktasks_mapper``::
else if (t->type == task_type_cooling || t->type == task_type_sourceterms) {
if (cell_is_active_hydro(t->ci, e)) scheduler_activate(s, t);
}
Then you will need to update the estimate for the number of tasks in ``engine_estimate_nr_tasks`` by modifying ``n1`` or ``n2``.
Initially, the engine will need to skip the task that updates the particles.
It is the case for the cooling, therefore you will need to add it in ``engine_skip_force_and_kick``.
Implementing your Task
----------------------
The last part is situated in ``runner.c``.
You will need to implement a function ``runner_do_cooling``
(do not forget to time it)::
void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
/* Now you can check if something is required at this time step.
* You may want to use a different cell_is_active function depending
* on your task
*/
if (!cell_is_active_hydro(c, e)) return;
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_cooling(r, c->progeny[k], 0);
} else {
/* Implement your cooling here */
}
if (timer) TIMER_TOC(timer_do_cooling);
}
and add a call to this function in ``runner_main``
in the switch::
case task_type_cooling:
runner_do_cooling(r, t->ci, 1);
break;
Finalizing your Task
--------------------
Now that you have done the easiest part, you can start debugging by implementing a test and/or an example.
Before creating your merge request with your new task, do not forget the most funny part that consists in writing a nice and beautiful documentation ;)
.. Task
Loic Hausammann 17th July 2018
.. _task:
Task System
===========
This section of the documentation includes information on the task system
available in SWIFT, as well as how to implement your own task.
SWIFT produces at the beginning of each simulation a ``dot`` file (see the graphviz library for more information).
It contains the full hierarchy of tasks used in this simulation.
You can convert the ``dot`` file into a ``png`` with the following command
``dot -Tpng dependency_graph.dot -o dependency_graph.png``.
.. toctree::
:maxdepth: 2
:caption: Contents:
adding_your_own
......@@ -20,3 +20,4 @@ difference is the parameter file that will need to be adapted for SWIFT.
Cooling/index
EquationOfState/index
NewOption/index
Task/index
......@@ -54,10 +54,12 @@ infile = args.input
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
"end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in",
"grav_down", "grav_mesh", "cooling", "sourceterms", "count"]
"grav_down", "grav_mesh", "cooling", "sourceterms",
"stars_ghost_in", "stars_ghost", "stars_ghost_out",
"count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
"tend", "xv", "rho", "gpart", "multipole", "spart", "stars_density", "count"]
SIDS = ["(-1,-1,-1)", "(-1,-1, 0)", "(-1,-1, 1)", "(-1, 0,-1)",
"(-1, 0, 0)", "(-1, 0, 1)", "(-1, 1,-1)", "(-1, 1, 0)",
......
......@@ -98,6 +98,7 @@ void print_help_message(void) {
printf(" %2s %14s %s\n", "-r", "", "Continue using restart files.");
printf(" %2s %14s %s\n", "-s", "", "Run with hydrodynamics.");
printf(" %2s %14s %s\n", "-S", "", "Run with stars.");
printf(" %2s %14s %s\n", "-b", "", "Run with stars feedback.");
printf(" %2s %14s %s\n", "-t", "{int}",
"The number of threads to use on each MPI rank. Defaults to 1 if not "
"specified.");
......@@ -135,6 +136,7 @@ int main(int argc, char *argv[]) {
struct gpart *gparts = NULL;
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
struct stars_props stars_properties;
struct part *parts = NULL;
struct phys_const prog_const;
struct sourceterms sourceterms;
......@@ -192,6 +194,7 @@ int main(int argc, char *argv[]) {
int with_self_gravity = 0;
int with_hydro = 0;
int with_stars = 0;
int with_feedback = 0;
int with_fp_exceptions = 0;
int with_drift_all = 0;
int with_mpole_reconstruction = 0;
......@@ -208,7 +211,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
while ((c = getopt(argc, argv, "acCdDef:FgGhMn:o:P:rsSt:Tv:xy:Y:")) != -1)
while ((c = getopt(argc, argv, "abcCdDef:FgGhMn:o:P:rsSt:Tv:xy:Y:")) != -1)
switch (c) {
case 'a':
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA)
......@@ -217,6 +220,9 @@ int main(int argc, char *argv[]) {
error("Need NUMA support for thread affinity");
#endif
break;
case 'b':
with_feedback = 1;
break;
case 'c':
with_cosmology = 1;
break;
......@@ -384,6 +390,15 @@ int main(int argc, char *argv[]) {
return 1;
}
if (!with_stars && with_feedback) {
if (myrank == 0)
printf(
"Error: Cannot process feedback without stars, -S must be "
"chosen.\n");
if (myrank == 0) print_help_message();
return 1;
}
/* Let's pin the main thread, now we know if affinity will be used. */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
if (with_aff &&
......@@ -681,6 +696,13 @@ int main(int argc, char *argv[]) {
else
bzero(&eos, sizeof(struct eos_parameters));
/* Initialise the stars properties */
if (with_stars)
stars_props_init(&stars_properties, &prog_const, &us, params,
&hydro_properties);
else
bzero(&stars_properties, sizeof(struct stars_props));
/* Initialise the gravity properties */
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology);
......@@ -752,7 +774,7 @@ int main(int argc, char *argv[]) {
/* Check once and for all that we don't have unwanted links */
if (!with_stars && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_star) error("Linking problem");
if (gparts[k].type == swift_type_stars) error("Linking problem");
}
if (!with_hydro && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
......@@ -778,7 +800,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0)
message(
"Read %lld gas particles, %lld star particles and %lld gparts from "
"Read %lld gas particles, %lld stars particles and %lld gparts from "
"the ICs.",
N_total[0], N_total[2], N_total[1]);
......@@ -887,6 +909,7 @@ int main(int argc, char *argv[]) {
if (with_cooling) engine_policies |= engine_policy_cooling;
if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
if (with_stars) engine_policies |= engine_policy_stars;
if (with_feedback) engine_policies |= engine_policy_feedback;
if (with_structure_finding)
engine_policies |= engine_policy_structure_finding;
......@@ -894,8 +917,8 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &gravity_properties, &mesh, &potential,
&cooling_func, &chemistry, &sourceterms);
&hydro_properties, &gravity_properties, &stars_properties,
&mesh, &potential, &cooling_func, &chemistry, &sourceterms);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
......@@ -910,7 +933,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
long long N_DM = N_total[1] - N_total[2] - N_total[0];
message(
"Running on %lld gas particles, %lld star particles and %lld DM "
"Running on %lld gas particles, %lld stars particles and %lld DM "
"particles (%lld gravity particles)",
N_total[0], N_total[2], N_total[1] > 0 ? N_DM : 0, N_total[1]);
message(
......
......@@ -55,6 +55,8 @@ Scheduler:
cell_sub_size_self_hydro: 32000 # (Optional) Maximal number of interactions per sub-self hydro task (this is the default value).
cell_sub_size_pair_grav: 256000000 # (Optional) Maximal number of interactions per sub-pair gravity task (this is the default value).
cell_sub_size_self_grav: 32000 # (Optional) Maximal number of interactions per sub-self gravity task (this is the default value).
cell_sub_size_pair_stars: 256000000 # (Optional) Maximal number of interactions per sub-pair stars task (this is the default value).
cell_sub_size_self_stars: 32000 # (Optional) Maximal number of interactions per sub-self stars task (this is the default value).
cell_split_size: 400 # (Optional) Maximal number of particles per cell (this is the default value).
cell_subdepth_grav: 2 # (Optional) Maximal depth the gravity tasks can be pushed down (this is the default value).
max_top_level_cells: 12 # (Optional) Maximal number of top-level cells in any dimension. The number of top-level cells will be the cube of this (this is the default value).
......
......@@ -112,17 +112,20 @@ pl.rcParams.update(PLOT_PARAMS)
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "init_grav_out", "ghost_in", "ghost", "ghost_out", "extra_ghost", "drift_part", "drift_gpart",
"end_force", "kick1", "kick2", "timestep", "send", "recv", "grav_long_range", "grav_mm", "grav_down_in",
"grav_down", "grav_mesh", "cooling", "sourceterms", "count"]
"grav_down", "grav_mesh", "cooling", "sourceterms",
"stars_ghost_in", "stars_ghost", "stars_ghost_out",
"count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
"tend", "xv", "rho", "gpart", "multipole", "spart", "stars_density", "count"]
# Task/subtypes of interest.
FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
"sub_self/density", "pair/force", "pair/density", "pair/grav",
"sub_pair/force",
"sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
"recv/tend", "send/tend", "recv/gpart", "send/gpart"]
"recv/tend", "send/tend", "recv/gpart", "send/gpart", "self/stars_density",
"pair/stars_density", "sub_self/stars_density", "sub_pair/stars_density"]
# A number of colours for the various types. Recycled when there are
# more task types than colours...
......
......@@ -126,8 +126,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
riemann/riemann_exact.h riemann/riemann_vacuum.h \
riemann/riemann_checks.h \
stars.h stars_io.h \
stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
stars/Default/star_debug.h stars/Default/star_part.h \
stars/Default/stars.h stars/Default/stars_iact.h stars/Default/stars_io.h \
stars/Default/stars_debug.h stars/Default/stars_part.h \
potential/none/potential.h potential/point_mass/potential.h \
potential/isothermal/potential.h potential/disc_patch/potential.h \
potential/sine_wave/potential.h \
......
......@@ -174,6 +174,21 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active_gravity(
return (c->ti_gravity_end_max == e->ti_current);
}
/**
* @brief Does a cell contain any s-particle finishing their time-step now ?
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell contains at least an active particle, 0 otherwise.