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

Commit 478c0b67 authored by Loic Hausammann's avatar Loic Hausammann Committed by Matthieu Schaller
Browse files

Sink cut off radius and first sink tasks

parent 11e3ab0e
......@@ -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;
}
}
}