gitlab is upgraded to version 13, please report any issues and enjoy

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

Merge branch 'sink_cut_off_radius' into 'master'

Sink cut off radius and first sink tasks

See merge request !1149
parents 11e3ab0e 478c0b67
......@@ -80,3 +80,6 @@ HernquistPotential:
bulgefraction: 0.014 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.2 # Softening size (internal units)
DefaultSink:
cut_off_radius: 0.1 # Cut off radius of all the sinks in internal units.
......@@ -96,6 +96,7 @@ int main(int argc, char *argv[]) {
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
struct stars_props stars_properties;
struct sink_props sink_properties;
struct feedback_props feedback_properties;
struct entropy_floor_properties entropy_floor;
struct black_holes_props black_holes_properties;
......@@ -988,6 +989,12 @@ int main(int argc, char *argv[]) {
} else
bzero(&black_holes_properties, sizeof(struct black_holes_props));
/* Initialise the sink properties */
if (with_sink) {
sink_props_init(&sink_properties, &prog_const, &us, params, &cosmo);
} else
bzero(&sink_properties, sizeof(struct sink_props));
/* Initialise the cooling function properties */
#ifdef COOLING_NONE
if (with_cooling) {
......@@ -1309,7 +1316,7 @@ int main(int argc, char *argv[]) {
N_total[swift_type_dark_matter_background], engine_policies,
talking, &reparttype, &us, &prog_const, &cosmo,
&hydro_properties, &entropy_floor, &gravity_properties,
&stars_properties, &black_holes_properties,
&stars_properties, &black_holes_properties, &sink_properties,
&feedback_properties, &mesh, &potential, &cooling_func,
&starform, &chemistry, &fof_properties, &los_properties);
engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
......
......@@ -561,3 +561,9 @@ EAGLEAGN:
threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
# Parameters related to the sink particles ---------------------------------------
# Default sink particles
DefaultSink:
cut_off_radius: 1e-3 # Cut off radius of the sink particles (in internal units). This parameter should be adapted with the resolution.
......@@ -295,8 +295,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h \
pressure_floor/GEAR/pressure_floor_iact.h pressure_floor/none/pressure_floor_iact.h \
pressure_floor/GEAR/pressure_floor_struct.h pressure_floor/none/pressure_floor_struct.h \
sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h \
sink.h sink_io.h
sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h sink/Default/sink_properties.h \
sink.h sink_io.h sink_properties.h
# Sources and special flags for the gravity library
......
......@@ -3608,11 +3608,13 @@ void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
/* Store the current dx_max and h_max values. */
ci->sinks.dx_max_part_old = ci->sinks.dx_max_part;
ci->sinks.r_cut_max_old = ci->sinks.r_cut_max;
ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
ci->hydro.h_max_old = ci->hydro.h_max;
if (cj != NULL) {
cj->sinks.dx_max_part_old = cj->sinks.dx_max_part;
cj->sinks.r_cut_max_old = cj->sinks.r_cut_max;
cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
cj->hydro.h_max_old = cj->hydro.h_max;
}
......@@ -4804,6 +4806,8 @@ int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
/* Unskip all the other task types. */
if (c->nodeID == nodeID && cell_is_active_sinks(c, e)) {
if (c->sinks.sink_in != NULL) scheduler_activate(s, c->sinks.sink_in);
if (c->sinks.sink_out != NULL) scheduler_activate(s, c->sinks.sink_out);
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
......@@ -5662,6 +5666,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
struct sink *const sinks = c->sinks.parts;
float dx_max = 0.f, dx2_max = 0.f;
float cell_r_max = 0.f;
/* Drift irrespective of cell flags? */
force = (force || cell_get_flag(c, cell_flag_do_sink_drift));
......@@ -5701,10 +5706,12 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
/* Update */
dx_max = max(dx_max, cp->sinks.dx_max_part);
cell_r_max = max(cell_r_max, cp->sinks.r_cut_max);
}
}
/* Store the values */
c->sinks.r_cut_max = cell_r_max;
c->sinks.dx_max_part = dx_max;
/* Update the time of the last drift */
......@@ -5769,12 +5776,17 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
}
}
/* sp->h does not need to be limited. */
/* Compute (square of) motion since last cell construction */
const float dx2 = sink->x_diff[0] * sink->x_diff[0] +
sink->x_diff[1] * sink->x_diff[1] +
sink->x_diff[2] * sink->x_diff[2];
dx2_max = max(dx2_max, dx2);
/* Maximal smoothing length */
cell_r_max = max(cell_r_max, sink->r_cut);
/* Get ready for a density calculation */
if (sink_is_active(sink, e)) {
sink_init_sink(sink);
......@@ -5785,6 +5797,7 @@ void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
dx_max = sqrtf(dx2_max);
/* Store the values */
c->sinks.r_cut_max = cell_r_max;
c->sinks.dx_max_part = dx_max;
/* Update the time of the last drift */
......
......@@ -766,6 +766,12 @@ struct cell {
/*! Nr of #sink this cell can hold after addition of new one. */
int count_total;
/*! Max cut off radius in this cell. */
float r_cut_max;
/*! Values of r_cut_max before the drifts, used for sub-cell tasks. */
float r_cut_max_old;
/*! Number of #sink updated in this cell. */
int updated;
......@@ -798,6 +804,13 @@ struct cell {
/*! The drift task for sinks */
struct task *drift;
/*! Implicit tasks marking the entry of the sink block of tasks
*/
struct task *sink_in;
/*! Implicit tasks marking the exit of the sink block of tasks */
struct task *sink_out;
} sinks;
#ifdef WITH_MPI
......@@ -1214,9 +1227,7 @@ __attribute__((always_inline)) INLINE static int
cell_can_recurse_in_self_sinks_task(const struct cell *c) {
/* Is the cell split and not smaller than the smoothing length? */
// loic TODO: add cut off radius
return c->split &&
//(kernel_gamma * c->sinks.h_max_old < 0.5f * c->dmin) &&
return c->split && (c->sinks.r_cut_max_old < 0.5f * c->dmin) &&
(kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
}
......
......@@ -191,6 +191,14 @@ int checkSpacehmax(struct space *s) {
}
}
float cell_sinks_h_max = 0.0f;
for (int k = 0; k < s->nr_cells; k++) {
if (s->cells_top[k].nodeID == s->e->nodeID &&
s->cells_top[k].sinks.r_cut_max > cell_sinks_h_max) {
cell_sinks_h_max = s->cells_top[k].sinks.r_cut_max;
}
}
/* Now all particles. */
float part_h_max = 0.0f;
for (size_t k = 0; k < s->nr_parts; k++) {
......@@ -207,9 +215,18 @@ int checkSpacehmax(struct space *s) {
}
}
/* Now all the sinks. */
float sink_h_max = 0.0f;
for (size_t k = 0; k < s->nr_sinks; k++) {
if (s->sinks[k].r_cut > sink_h_max) {
sink_h_max = s->sinks[k].r_cut;
}
}
/* If within some epsilon we are OK. */
if (fabsf(cell_h_max - part_h_max) <= FLT_EPSILON &&
fabsf(cell_stars_h_max - spart_h_max) <= FLT_EPSILON)
fabsf(cell_stars_h_max - spart_h_max) <= FLT_EPSILON &&
fabsf(cell_sinks_h_max - sink_h_max) <= FLT_EPSILON)
return 1;
/* There is a problem. Hunt it down. */
......@@ -247,6 +264,23 @@ int checkSpacehmax(struct space *s) {
}
}
/* sink */
for (int k = 0; k < s->nr_cells; k++) {
if (s->cells_top[k].nodeID == s->e->nodeID) {
if (s->cells_top[k].sinks.r_cut_max > sink_h_max) {
message("cell %d is inconsistent (%f > %f)", k,
s->cells_top[k].sinks.r_cut_max, sink_h_max);
}
}
}
for (size_t k = 0; k < s->nr_sinks; k++) {
if (s->sinks[k].r_cut > cell_sinks_h_max) {
message("spart %lld is inconsistent (%f > %f)", s->sinks[k].id,
s->sinks[k].r_cut, cell_sinks_h_max);
}
}
return 0;
}
......@@ -267,6 +301,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
float dx_max = 0.0f;
float stars_h_max = 0.0f;
float stars_dx_max = 0.0f;
float sinks_h_max = 0.0f;
float sinks_dx_max = 0.0f;
int result = 1;
const double loc_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
......@@ -329,6 +365,33 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
stars_dx_max = max(stars_dx_max, sqrt(dx2));
}
const size_t nr_sinks = c->sinks.count;
struct sink *sinks = c->sinks.parts;
for (size_t k = 0; k < nr_sinks; k++) {
struct sink *const sp = &sinks[k];
if (sp->x[0] < loc_min[0] || sp->x[0] >= loc_max[0] ||
sp->x[1] < loc_min[1] || sp->x[1] >= loc_max[1] ||
sp->x[2] < loc_min[2] || sp->x[2] >= loc_max[2]) {
message(
"Inconsistent sink position p->x=[%e %e %e], c->loc=[%e %e %e] "
"c->width=[%e %e %e]",
sp->x[0], sp->x[1], sp->x[2], c->loc[0], c->loc[1], c->loc[2],
c->width[0], c->width[1], c->width[2]);
result = 0;
}
const float dx2 = sp->x_diff[0] * sp->x_diff[0] +
sp->x_diff[1] * sp->x_diff[1] +
sp->x_diff[2] * sp->x_diff[2];
sinks_h_max = max(sinks_h_max, sp->r_cut);
sinks_dx_max = max(sinks_dx_max, sqrt(dx2));
}
if (c->split) {
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL) {
......@@ -365,6 +428,19 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
result = 0;
}
if (c->sinks.r_cut_max != sinks_h_max) {
message("%d Inconsistent sinks_h_max: cell %f != parts %f", *depth,
c->sinks.r_cut_max, sinks_h_max);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
if (c->sinks.dx_max_part != sinks_dx_max) {
message("%d Inconsistent stars_dx_max: %f != %f", *depth,
c->sinks.dx_max_part, sinks_dx_max);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
return result;
}
......
......@@ -89,6 +89,7 @@
#include "runner.h"
#include "serial_io.h"
#include "single_io.h"
#include "sink_properties.h"
#include "sort_part.h"
#include "star_formation.h"
#include "star_formation_logger.h"
......@@ -2078,7 +2079,8 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_stars_resort || t->type == task_type_extra_ghost ||
t->type == task_type_stars_ghost ||
t->type == task_type_stars_ghost_in ||
t->type == task_type_stars_ghost_out ||
t->type == task_type_stars_ghost_out || t->type == task_type_sink_in ||
t->type == task_type_sink_out ||
t->type == task_type_bh_swallow_ghost1 ||
t->type == task_type_bh_swallow_ghost2 ||
t->type == task_type_bh_swallow_ghost3 || t->type == task_type_bh_in ||
......@@ -2461,6 +2463,20 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
}
}
if (s->cells_top != NULL && s->nr_sinks > 0) {
for (int i = 0; i < s->nr_cells; i++) {
struct cell *c = &s->cells_top[i];
if (c->nodeID == engine_rank && c->sinks.count > 0) {
float sink_h_max = c->sinks.parts[0].r_cut;
for (int k = 1; k < c->sinks.count; k++) {
if (c->sinks.parts[k].r_cut > sink_h_max)
sink_h_max = c->sinks.parts[k].r_cut;
}
c->sinks.r_cut_max = max(sink_h_max, c->sinks.r_cut_max);
}
}
}
clocks_gettime(&time2);
#ifdef SWIFT_DEBUG_CHECKS
......@@ -3881,6 +3897,7 @@ static void engine_dumper_init(struct engine *e) {
* @param gravity The #gravity_props used for this run.
* @param stars The #stars_props used for this run.
* @param black_holes The #black_holes_props used for this run.
* @param sinks The #sink_props used for this run.
* @param feedback The #feedback_props used for this run.
* @param mesh The #pm_mesh used for the long-range periodic forces.
* @param potential The properties of the external potential.
......@@ -3901,6 +3918,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
const struct black_holes_props *black_holes,
const struct sink_props *sinks,
struct feedback_props *feedback, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
......@@ -3982,6 +4000,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->gravity_properties = gravity;
e->stars_properties = stars;
e->black_holes_properties = black_holes;
e->sink_properties = sinks;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
......@@ -5442,6 +5461,7 @@ void engine_clean(struct engine *e, const int fof, const int restart) {
free((void *)e->output_options);
free((void *)e->external_potential);
free((void *)e->black_holes_properties);
free((void *)e->sink_properties);
free((void *)e->stars_properties);
free((void *)e->gravity_properties);
free((void *)e->hydro_properties);
......@@ -5509,6 +5529,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
starformation_struct_dump(e->star_formation, stream);
feedback_struct_dump(e->feedback_props, stream);
black_holes_struct_dump(e->black_holes_properties, stream);
sink_struct_dump(e->sink_properties, stream);
chemistry_struct_dump(e->chemistry, stream);
#ifdef WITH_FOF
fof_struct_dump(e->fof_properties, stream);
......@@ -5634,6 +5655,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
black_holes_struct_restore(black_holes_properties, stream);
e->black_holes_properties = black_holes_properties;
struct sink_props *sink_properties =
(struct sink_props *)malloc(sizeof(struct sink_props));
sink_struct_restore(sink_properties, stream);
e->sink_properties = sink_properties;
struct chemistry_global_data *chemistry =
(struct chemistry_global_data *)malloc(
sizeof(struct chemistry_global_data));
......
......@@ -447,6 +447,9 @@ struct engine {
/* Properties of the black hole model */
const struct black_holes_props *black_holes_properties;
/* Properties of the sink model */
const struct sink_props *sink_properties;
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
......@@ -569,6 +572,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct entropy_floor_properties *entropy_floor,
struct gravity_props *gravity, const struct stars_props *stars,
const struct black_holes_props *black_holes,
const struct sink_props *sinks,
struct feedback_props *feedback, struct pm_mesh *mesh,
const struct external_potential *potential,
struct cooling_function_data *cooling_func,
......
......@@ -1239,6 +1239,25 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
}
}
/* Subgrid tasks: sink */
if (with_sinks) {
c->sinks.sink_in =
scheduler_addtask(s, task_type_sink_in, task_subtype_none, 0,
/* implicit = */ 1, c, NULL);
c->sinks.sink_out =
scheduler_addtask(s, task_type_sink_out, task_subtype_none, 0,
/* implicit = */ 1, c, NULL);
/* TODO */
scheduler_addunlock(s, c->sinks.sink_in, c->sinks.sink_out);
/* Link to the main tasks */
scheduler_addunlock(s, c->super->kick2, c->sinks.sink_in);
scheduler_addunlock(s, c->sinks.sink_out, c->super->timestep);
}
/* Subgrid tasks: black hole feedback */
if (with_black_holes) {
......
......@@ -947,6 +947,12 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate(s, t);
}
/* Feedback implicit tasks? */
else if (t_type == task_type_sink_in || t_type == task_type_sink_out) {
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 ||
......
......@@ -24,6 +24,7 @@
/* Local includes */
#include "minmax.h"
#include "sink_part.h"
#include "sink_properties.h"
/**
* @brief Computes the time-step of a given sink particle.
......@@ -45,8 +46,9 @@ __attribute__((always_inline)) INLINE static float sink_compute_timestep(
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void sink_first_init_sink(
struct sink* sp) {
struct sink* sp, const struct sink_props* sink_props) {
sp->r_cut = sink_props->cut_off_radius;
sp->time_bin = 0;
}
......
......@@ -43,6 +43,9 @@ struct sink {
/*! Particle velocity. */
float v[3];
/*! Cut off radius. */
float r_cut;
/*! Sink particle mass */
float mass;
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
*
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_DEFAULT_SINK_PROPERTIES_H
#define SWIFT_DEFAULT_SINK_PROPERTIES_H
/**
* @brief Properties of sink in the Default model.
*/
struct sink_props {
/*! Cut off radius */
float cut_off_radius;
};
/**
* @brief Initialise the sink properties from the parameter file.
*
* @param sp The #sink_props.
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
* @param params The parsed parameters.
* @param cosmo The cosmological model.
*/
INLINE static void sink_props_init(struct sink_props *sp,
const struct phys_const *phys_const,
const struct unit_system *us,
struct swift_params *params,
const struct cosmology *cosmo) {
sp->cut_off_radius =
parser_get_param_float(params, "DefaultSink:cut_off_radius");
}
/**
* @brief Write a sink_props struct to the given FILE as a stream of
* bytes.
*
* @param props the sink properties struct
* @param stream the file stream
*/
INLINE static void sink_struct_dump(const struct sink_props *props,
FILE *stream) {
restart_write_blocks((void *)props, sizeof(struct sink_props), 1, stream,
"sink props", "Sink props");
}
/**
* @brief Restore a sink_props struct from the given FILE as a stream of
* bytes.
*
* @param props the sink properties struct
* @param stream the file stream
*/
INLINE static void sink_struct_restore(const struct sink_props *props,
FILE *stream) {
restart_read_blocks((void *)props, sizeof(struct sink_props), 1, stream, NULL,
"Sink props");
}
#endif /* SWIFT_DEFAULT_SINK_PROPERTIES_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
*
* 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/>.
*
******************************************************************************/
#ifndef SWIFT_SINK_PROPERTIES_H
#define SWIFT_SINK_PROPERTIES_H
/* Config parameters. */
#include "../config.h"
/* Select the correct sink model */
#if defined(SINK_NONE)
#include "./sink/Default/sink_properties.h"
#else
#error "Invalid choice of black hole model"
#endif
#endif /* SWIFT_SINK_PROPERTIES_H */
......@@ -260,6 +260,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->black_holes.drift = NULL;
c->black_holes.black_holes_in = NULL;
c->black_holes.black_holes_out = NULL;
c->sinks.sink_in = NULL;
c->sinks.sink_out = NULL;
c->grav.drift = NULL;
c->grav.drift_out = NULL;
c->hydro.cooling_in = NULL;
......@@ -380,6 +382,7 @@ void space_regrid(struct space *s, int verbose) {
const size_t nr_parts = s->nr_parts;
const size_t nr_sparts = s->nr_sparts;
const size_t nr_bparts = s->nr_bparts;
const size_t nr_sinks = s->nr_sinks;
const ticks tic = getticks();
const integertime_t ti_current = (s->e != NULL) ? s->e->ti_current : 0;
......@@ -402,6 +405,9 @@ void space_regrid(struct space *s, int verbose) {
if (c->black_holes.h_max > h_max) {
h_max = c->black_holes.h_max;
}
if (c->sinks.r_cut_max > h_max) {
h_max = c->sinks.r_cut_max / kernel_gamma;
}
}
/* Can we instead use all the top-level cells? */
......@@ -417,6 +423,9 @@ void space_regrid(struct space *s, int verbose) {
if (c->nodeID == engine_rank && c->black_holes.h_max > h_max) {
h_max = c->black_holes.h_max;
}
if (c->nodeID == engine_rank && c->sinks.r_cut_max > h_max) {
h_max = c->sinks.r_cut_max / kernel_gamma;
}
}
/* Last option: run through the particles */
......@@ -430,6 +439,9 @@ void space_regrid(struct space *s, int verbose) {
for (size_t k = 0; k < nr_bparts; k++) {
if (s->bparts[k].h > h_max) h_max = s->bparts[k].h;
}
for (size_t k = 0; k < nr_sinks; k++) {
if (s->sinks[k].r_cut > h_max) h_max = s->sinks[k].r_cut / kernel_gamma;
}
}
}