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

Compute the sftening lengths from the DM particle masses that have been read in.

parent 37de2088
......@@ -745,12 +745,6 @@ int main(int argc, char *argv[]) {
const int generate_gas_in_ics = parser_get_opt_param_int(
params, "InitialConditions:generate_gas_in_ics", 0);
/* Some checks that we are not doing something stupid */
if (generate_gas_in_ics && flag_entropy_ICs)
error("Can't generate gas if the entropy flag is set in the ICs.");
if (generate_gas_in_ics && !with_cosmology)
error("Can't generate gas if the run is not cosmological.");
/* Initialise the cosmology */
if (with_cosmology)
cosmology_init(params, &us, &prog_const, &cosmo);
......@@ -804,12 +798,6 @@ int main(int argc, char *argv[]) {
} else
bzero(&black_holes_properties, sizeof(struct black_holes_props));
/* Initialise the gravity properties */
bzero(&gravity_properties, sizeof(struct gravity_props));
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
with_cosmology, periodic);
/* Initialise the cooling function properties */
#ifdef COOLING_NONE
if (with_cooling || with_temperature) {
......@@ -867,24 +855,26 @@ int main(int argc, char *argv[]) {
#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run);
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, &high_res_DM_mass, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
nr_threads, dry_run);
#else
read_ic_serial(
ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts, &Ngas, &Ngpart,
&Ngpart_background, &Nspart, &Nbpart, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h, cosmo.a, myrank,
nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads, dry_run);
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, &high_res_DM_mass, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
cosmo.a, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
nr_threads, dry_run);
#endif
#else
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &bparts,
&Ngas, &Ngpart, &Ngpart_background, &Nspart, &Nbpart,
&flag_entropy_ICs, with_hydro,
&flag_entropy_ICs, &high_res_DM_mass, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
with_black_holes, cleanup_h, cleanup_sqrt_a, cosmo.h,
cosmo.a, nr_threads, dry_run);
......@@ -897,6 +887,12 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
/* Some checks that we are not doing something stupid */
if (generate_gas_in_ics && flag_entropy_ICs)
error("Can't generate gas if the entropy flag is set in the ICs.");
if (generate_gas_in_ics && !with_cosmology)
error("Can't generate gas if the run is not cosmological.");
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
if (!with_stars && !dry_run) {
......@@ -970,6 +966,12 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
/* Initialise the gravity properties */
bzero(&gravity_properties, sizeof(struct gravity_props));
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, high_res_DM_mass,
with_cosmology, periodic);
/* Initialise the external potential properties */
bzero(&potential, sizeof(struct external_potential));
if (with_external_gravity)
......
......@@ -2175,8 +2175,11 @@ void cell_reset_task_counters(struct cell *c) {
*
* @param c The #cell.
* @param ti_current The current integer time.
* @param grav_props The properties of the gravity scheme.
*/
void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
void cell_make_multipoles(struct cell *c, integertime_t ti_current,
const struct gravity_props *const grav_props) {
/* Reset everything */
gravity_reset(c->grav.multipole);
......@@ -2184,7 +2187,7 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
/* Start by recursing */
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL)
cell_make_multipoles(c->progeny[k], ti_current);
cell_make_multipoles(c->progeny[k], ti_current, grav_props);
}
/* Compute CoM of all progenies */
......@@ -2245,7 +2248,7 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
} else {
if (c->grav.count > 0) {
gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count);
gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count, grav_props);
const double dx =
c->grav.multipole->CoM[0] > c->loc[0] + c->width[0] * 0.5
? c->grav.multipole->CoM[0] - c->loc[0]
......@@ -2322,7 +2325,9 @@ void cell_check_foreign_multipole(const struct cell *c) {
*
* @param c Cell to act upon
*/
void cell_check_multipole(struct cell *c) {
void cell_check_multipole(struct cell *c,
const struct gravity_props *const grav_props) {
#ifdef SWIFT_DEBUG_CHECKS
struct gravity_tensors ma;
const double tolerance = 1e-3; /* Relative */
......@@ -2330,11 +2335,12 @@ void cell_check_multipole(struct cell *c) {
/* First recurse */
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k]);
if (c->progeny[k] != NULL)
cell_check_multipole(c->progeny[k], grav_props);
if (c->grav.count > 0) {
/* Brute-force calculation */
gravity_P2M(&ma, c->grav.parts, c->grav.count);
gravity_P2M(&ma, c->grav.parts, c->grav.count, grav_props);
/* Now compare the multipole expansion */
if (!gravity_multipole_equal(&ma, c->grav.multipole, tolerance)) {
......
......@@ -852,8 +852,10 @@ int cell_link_foreign_gparts(struct cell *c, struct gpart *gparts);
int cell_count_parts_for_tasks(const struct cell *c);
int cell_count_gparts_for_tasks(const struct cell *c);
void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c);
void cell_make_multipoles(struct cell *c, integertime_t ti_current,
const struct gravity_props *const grav_props);
void cell_check_multipole(struct cell *c,
const struct gravity_props *const grav_props);
void cell_check_foreign_multipole(const struct cell *c);
void cell_clean(struct cell *c);
void cell_check_part_drift_point(struct cell *c, void *data);
......
......@@ -4198,7 +4198,7 @@ void engine_do_reconstruct_multipoles_mapper(void *map_data, int num_elements,
if (c != NULL && c->nodeID == e->nodeID) {
/* Construct the multipoles in this cell hierarchy */
cell_make_multipoles(c, e->ti_current);
cell_make_multipoles(c, e->ti_current, e->gravity_properties);
}
}
}
......
......@@ -48,7 +48,10 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
__attribute__((always_inline)) INLINE static float gravity_get_softening(
const struct gpart* gp, const struct gravity_props* restrict grav_props) {
return grav_props->epsilon_cur;
if (gp->type == swift_type_dark_matter_background)
return grav_props->epsilon_background_fac * cbrtf(gp->mass);
else
return grav_props->epsilon_cur;
}
/**
......
......@@ -67,6 +67,11 @@ INLINE static void convert_gpart_vel(const struct engine* e,
ret[2] *= cosmo->a_inv;
}
INLINE static void convert_gpart_soft(const struct engine* e,
const struct gpart* gp, float* ret) {
ret[0] = gravity_get_softening(gp, e->gravity_properties);
}
/**
* @brief Specifies which g-particle fields to read from a dataset
*
......@@ -104,7 +109,7 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
int* num_fields) {
/* Say how much we want to write */
*num_fields = 5;
*num_fields = 6;
/* List what we want to write */
list[0] = io_make_output_field_convert_gpart(
......@@ -117,6 +122,8 @@ INLINE static void darkmatter_write_particles(const struct gpart* gparts,
UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
list[4] = io_make_output_field("GroupIDs", INT, 1, UNIT_CONV_NO_UNITS, gparts,
group_id);
list[5] = io_make_output_field_convert_gpart(
"Softenings", FLOAT, 1, UNIT_CONV_LENGTH, gparts, convert_gpart_soft);
}
#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
......@@ -44,8 +44,7 @@ struct gpart {
/*! Type of the #gpart (DM, gas, star, ...) */
enum part_type type;
/* Particle group ID and size in the FOF. */
size_t group_id, group_size;
unsigned char group_id, group_size;
#ifdef SWIFT_DEBUG_CHECKS
......
......@@ -40,8 +40,10 @@
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
const struct cosmology *cosmo, int with_cosmology,
int periodic) {
const struct cosmology *cosmo,
const double high_res_DM_mass,
const int with_cosmology,
const int periodic) {
/* Tree updates */
p->rebuild_frequency =
......@@ -89,10 +91,39 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
/* Softening parameters */
if (with_cosmology) {
p->epsilon_comoving =
parser_get_param_double(params, "Gravity:comoving_softening");
/* Check that a mass for the high res. particles was indeed read. */
if (high_res_DM_mass == -1.)
error("DM particle mass for type 1 has not been read from the ICs.");
const double ratio =
parser_get_param_double(params, "Gravity:softening_ratio");
/* Compute the comoving softening length as a fraction of the mean
* inter-particle density of the high-res. DM particles
* and assuming the mean density of the Universe is used in the simulation.
*/
const double mean_matter_density =
cosmo->Omega_m * cosmo->critical_density_0;
p->epsilon_comoving = ratio * cbrt(high_res_DM_mass / mean_matter_density);
/* Maximal physical softening taken straight from the parameter file */
p->epsilon_max_physical =
parser_get_param_double(params, "Gravity:max_physical_softening");
/* Compute the comoving softening length for background particles.
* Since they have variable masses the mass factor will be multiplied
* in later on.
* Note that we already multiply in the conversion from Plummer -> real
* softening length */
const double ratio_background =
parser_get_param_double(params, "Gravity:softening_ratio_background");
p->epsilon_background_fac = kernel_gravity_softening_plummer_equivalent *
ratio_background *
cbrt(1. / mean_matter_density);
} else {
p->epsilon_max_physical =
parser_get_param_double(params, "Gravity:max_physical_softening");
......@@ -109,7 +140,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
void gravity_props_update(struct gravity_props *p,
const struct cosmology *cosmo) {
/* Current softening lengths */
/* Current softening length for the high-resolution particles. */
double softening;
if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
softening = p->epsilon_max_physical / cosmo->a;
......@@ -121,11 +152,6 @@ void gravity_props_update(struct gravity_props *p,
/* Store things */
p->epsilon_cur = softening;
/* Other factors */
p->epsilon_cur2 = softening * softening;
p->epsilon_cur_inv = 1. / softening;
p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
}
void gravity_props_print(const struct gravity_props *p) {
......@@ -143,14 +169,14 @@ void gravity_props_print(const struct gravity_props *p) {
kernel_gravity_softening_name);
message(
"Self-gravity comoving softening: epsilon=%.4f (Plummer equivalent: "
"%.4f)",
"Self-gravity comoving softening: epsilon=%.6f (Plummer equivalent: "
"%.6f)",
p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
p->epsilon_comoving);
message(
"Self-gravity maximal physical softening: epsilon=%.4f (Plummer "
"equivalent: %.4f)",
"Self-gravity maximal physical softening: epsilon=%.6f (Plummer "
"equivalent: %.6f)",
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
p->epsilon_max_physical);
......@@ -171,33 +197,38 @@ void gravity_props_print(const struct gravity_props *p) {
void gravity_props_print_snapshot(hid_t h_grpgrav,
const struct gravity_props *p) {
io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
io_write_attribute_s(h_grpgrav, "Softening style",
kernel_gravity_softening_name);
io_write_attribute_f(
h_grpgrav, "Comoving softening length [internal units]",
p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(
h_grpgrav,
"Comoving Softening length (Plummer equivalent) [internal units]",
p->epsilon_comoving);
io_write_attribute_f(
h_grpgrav, "Maximal physical softening length [internal units]",
p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
io_write_attribute_f(h_grpgrav,
"Maximal physical softening length (Plummer equivalent) "
" [internal units]",
p->epsilon_max_physical);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION);
io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio", p->r_cut_max_ratio);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio", p->r_cut_min_ratio);
io_write_attribute_f(h_grpgrav, "Tree update frequency",
p->rebuild_frequency);
io_write_attribute_s(h_grpgrav, "Mesh truncation function",
kernel_long_gravity_truncation_name);
/* io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta); */
/* io_write_attribute_s(h_grpgrav, "Softening style", */
/* kernel_gravity_softening_name); */
/* io_write_attribute_f( */
/* h_grpgrav, "Comoving softening length [internal units]", */
/* p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent); */
/* io_write_attribute_f( */
/* h_grpgrav, */
/* "Comoving Softening length (Plummer equivalent) [internal units]", */
/* p->epsilon_comoving); */
/* io_write_attribute_f( */
/* h_grpgrav, "Maximal physical softening length [internal units]", */
/* p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
*/
/* io_write_attribute_f(h_grpgrav, */
/* "Maximal physical softening length (Plummer
* equivalent) " */
/* " [internal units]", */
/* p->epsilon_max_physical); */
/* io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit); */
/* io_write_attribute_s(h_grpgrav, "Scheme", GRAVITY_IMPLEMENTATION); */
/* io_write_attribute_i(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
*/
/* io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth); */
/* io_write_attribute_f(h_grpgrav, "Mesh r_cut_max ratio",
* p->r_cut_max_ratio); */
/* io_write_attribute_f(h_grpgrav, "Mesh r_cut_min ratio",
* p->r_cut_min_ratio); */
/* io_write_attribute_f(h_grpgrav, "Tree update frequency", */
/* p->rebuild_frequency); */
/* io_write_attribute_s(h_grpgrav, "Mesh truncation function", */
/* kernel_long_gravity_truncation_name); */
}
#endif
......
......@@ -40,25 +40,18 @@ struct swift_params;
*/
struct gravity_props {
/*! Frequency of tree-rebuild in units of #gpart updates. */
float rebuild_frequency;
/* -------------- Softening for the regular DM and baryons ----------- */
/*! Periodic long-range mesh side-length */
int mesh_size;
/*! Softening length for high-res. particles at the current redshift. */
float epsilon_cur;
/*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
/* -------------- Softening for the background DM -------------------- */
/*! Distance below which the truncated mesh force is Newtonian in units of
* a_smooth */
float r_cut_min_ratio;
/*! Conversion factor from cbrt of particle mass to softening assuming
* a constant fraction of the mean inter-particle separation at that mass. */
float epsilon_background_fac;
/*! Distance above which the truncated mesh force is negligible in units of
* a_smooth */
float r_cut_max_ratio;
/*! Time integration dimensionless multiplier */
float eta;
/* -------------- Properties of the FFM gravity ---------------------- */
/*! Tree opening angle (Multipole acceptance criterion) */
double theta_crit;
......@@ -69,23 +62,44 @@ struct gravity_props {
/*! Inverse of opening angle */
double theta_crit_inv;
/*! Comoving softening */
double epsilon_comoving;
/* ------------- Properties of the softened gravity ------------------ */
/*! Maxium physical softening */
double epsilon_max_physical;
/*! Fraction of the mean inter particle separation corresponding to the
* co-moving softening length of the low-res. particles (DM + baryons) */
float mean_inter_particle_fraction_high_res;
/*! Current softening length */
float epsilon_cur;
/*! Co-moving softening length for for high-res. particles (DM + baryons)
* assuming a constant fraction of the mean inter-particle separation
* and a constant particle mass */
float epsilon_comoving;
/*! Square of current softening length */
float epsilon_cur2;
/*! Maximal softening length in physical coordinates for the high-res.
* particles (DM + baryons) */
float epsilon_max_physical;
/*! Inverse of current softening length */
float epsilon_cur_inv;
/* ------------- Properties of the time integration ----------------- */
/*! Cube of the inverse of current oftening length */
float epsilon_cur_inv3;
/*! Frequency of tree-rebuild in units of #gpart updates. */
float rebuild_frequency;
/*! Time integration dimensionless multiplier */
float eta;
/* ------------- Properties of the mesh-based gravity ---------------- */
/*! Periodic long-range mesh side-length */
int mesh_size;
/*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
/*! Distance below which the truncated mesh force is Newtonian in units of
* a_smooth */
float r_cut_min_ratio;
/*! Distance above which the truncated mesh force is negligible in units of
* a_smooth */
float r_cut_max_ratio;
/*! Gravitational constant (in internal units, copied from the physical
* constants) */
......@@ -94,9 +108,10 @@ struct gravity_props {
void gravity_props_print(const struct gravity_props *p);
void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct phys_const *phys_const,
const struct cosmology *cosmo, int with_cosmology,
int periodic);
const struct phys_const *phys_const,
const struct cosmology *cosmo,
const double high_res_DM_mass, const int with_cosmology,
const int periodic);
void gravity_props_update(struct gravity_props *p,
const struct cosmology *cosmo);
......
......@@ -338,7 +338,8 @@ void hydro_props_update(struct hydro_props *p, const struct gravity_props *gp,
* is a fixed fraction of the radius at which the softened forces
* recover a Newtonian behaviour (i.e. 2.8 * Plummer equivalent softening
* in the case of a cubic spline kernel). */
p->h_min = p->h_min_ratio * gp->epsilon_cur / kernel_gamma;
// MATTHIEU
p->h_min = p->h_min_ratio * 1. /*gp->epsilon_cur*/ / kernel_gamma;
}
/**
......
......@@ -117,6 +117,9 @@ struct multipole {
/*! Minimal velocity along each axis of all #gpart */
float min_delta_vel[3];
/*! Maximal co-moving softening of all the #gpart in the mulipole */
float max_softening;
/* 0th order term */
float M_000;
......@@ -459,6 +462,7 @@ INLINE static void gravity_multipole_init(struct multipole *m) {
*/
INLINE static void gravity_multipole_print(const struct multipole *m) {
printf("eps_max = %12.5e\n", m->max_softening);
printf("Vel= [%12.5e %12.5e %12.5e]\n", m->vel[0], m->vel[1], m->vel[2]);
printf("-------------------------\n");
printf("M_000= %12.5e\n", m->M_000);
......@@ -512,6 +516,9 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_add(struct multipole *restrict ma,
const struct multipole *restrict mb) {
/* Maximum of both softenings */
ma->max_softening = max(ma->max_softening, mb->max_softening);
/* Add 0th order term */
ma->M_000 += mb->M_000;
......@@ -630,6 +637,14 @@ INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
const double v2 = ma->vel[0] * ma->vel[0] + ma->vel[1] * ma->vel[1] +
ma->vel[2] * ma->vel[2];
/* Check maximal softening */
if (fabsf(ma->max_softening - mb->max_softening) /
fabsf(ma->max_softening + mb->max_softening) >
tolerance) {
message("max softening different!");
return 0;
}
/* Check bulk velocity (if non-zero and component > 1% of norm)*/
if (fabsf(ma->vel[0] + mb->vel[0]) > 1e-10 &&
(ma->vel[0] * ma->vel[0]) > 0.0001 * v2 &&
......@@ -1022,11 +1037,14 @@ INLINE static int gravity_multipole_equal(const struct gravity_tensors *ga,
* @param multi The #multipole (content will be overwritten).
* @param gparts The #gpart.
* @param gcount The number of particles.
* @param grav_props The properties of the gravity scheme.
*/
INLINE static void gravity_P2M(struct gravity_tensors *multi,
const struct gpart *gparts, int gcount) {
const struct gpart *gparts, const int gcount,
const struct gravity_props *const grav_props) {
/* Temporary variables */
float epsilon_max = 0.f;
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
double vel[3] = {0.f, 0.f, 0.f};
......@@ -1034,12 +1052,14 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
/* Collect the particle data for CoM. */
for (int k = 0; k < gcount; k++) {
const double m = gparts[k].mass;
const float epsilon = gravity_get_softening(&gparts[k], grav_props);
#ifdef SWIFT_DEBUG_CHECKS
if (gparts[k].time_bin == time_bin_inhibited)
error("Inhibited particle in P2M. Should have been removed earlier.");
#endif
epsilon_max = max(epsilon_max, epsilon);
mass += m;
com[0] += gparts[k].x[0] * m;
com[1] += gparts[k].x[1] * m;
......@@ -1203,6 +1223,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *multi,
#endif
/* Store the data on the multipole. */
multi->m_pole.max_softening = epsilon_max;
multi->m_pole.M_000 = mass;
multi->r_max = sqrt(r_max2);
multi->CoM[0] = com[0];
......@@ -1316,6 +1337,9 @@ INLINE static void gravity_M2M(struct multipole *restrict m_a,
const struct multipole *restrict m_b,
const double pos_a[3], const double pos_b[3]) {
/* "shift" the softening */
m_a->max_softening = m_b->max_softening;
/* Shift 0th order term */
m_a->M_000 = m_b->M_000;
......@@ -1964,8 +1988,8 @@ INLINE static void gravity_M2L_nonsym(
const int periodic, const double dim[3], const float rs_inv) {
/* Recover some constants */
const float eps = props->epsilon_cur;
const float eps_inv = props->epsilon_cur_inv;
const float eps = m_a->max_softening;
const float eps_inv = 1.f / eps;
/* Compute distance vector */
float dx = (float)(pos_b[0] - pos_a[0]);
......@@ -2015,8 +2039,8 @@ INLINE static void gravity_M2L_symmetric(
const float rs_inv) {
/* Recover some constants */
const float eps = props->epsilon_cur;
const float eps_inv = props->epsilon_cur_inv;
const float eps = m_a->max_softening;
const float eps_inv = 1.f / eps;
/* Compute distance vector */
float dx = (float)(pos_b[0] - pos_a[0]);
......
......@@ -493,11 +493,11 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
size_t* Ngparts, size_t* Ngparts_background, size_t* Nstars,
size_t* Nblackholes, int* flag_entropy, int with_hydro,
int with_gravity, int with_stars, int with_black_holes,
int cleanup_h, int cleanup_sqrt_a, double h, double a,
int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,