diff --git a/examples/main.c b/examples/main.c
index e795e210ca84b87410e1503c063d30c9a114ae83..56664f8dde5c6bb2667bfed6181e480b52977bd2 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -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)
diff --git a/src/cell.c b/src/cell.c
index 80c5798c2b5a13bcb3c283eb16fa58d537d63bf6..dfae80d5a2cf4a46fe8bc42052220be559efe31b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -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)) {
diff --git a/src/cell.h b/src/cell.h
index 9c9e2e64f85305e37fe912e0436e5dbde333f6d8..c6d32c8219b31cadfd20a7110226aa7e15badcf9 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -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);
diff --git a/src/engine.c b/src/engine.c
index 8f14cfd0088cff9d47bee38662a5aeff357dd8e0..4afbeb9b2f891d339d1c4d96d072e0b2c5a525ec 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -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);
}
}
}
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
index 6d1270c952100f3a25202fcdb22be09f9acaa8d9..aae8c9ee670019af3bf7f82cad0f75876810df30 100644
--- a/src/gravity/MultiSoftening/gravity.h
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -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;
}
/**
diff --git a/src/gravity/MultiSoftening/gravity_io.h b/src/gravity/MultiSoftening/gravity_io.h
index 2e443e8fdc2479f3f2feff30281dccad9a7b6519..11c970aa2973996916060bc38ffcccd8a3fe2fcc 100644
--- a/src/gravity/MultiSoftening/gravity_io.h
+++ b/src/gravity/MultiSoftening/gravity_io.h
@@ -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 */
diff --git a/src/gravity/MultiSoftening/gravity_part.h b/src/gravity/MultiSoftening/gravity_part.h
index 2c16a34a11d3a7674f5a6c122da780aeaff4a9dc..276c7f5176002b8a2e51768c47c6c91d8187fe51 100644
--- a/src/gravity/MultiSoftening/gravity_part.h
+++ b/src/gravity/MultiSoftening/gravity_part.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
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 8d3dfe9c4683842b905959ab3c688cc1f23309e9..c1807630dda831e45b63bfdea73a3a66750bed2f 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -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
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 79d72765e767351e29ffa30b615c1f7f4f98022c..81ec91ee1e0813b8779caccd540f87088ada8258 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -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);
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 0082546ce14030a37fa63d87ea88bc99153b8213..d823675579a5e13e5f80f2c5030a755d6e96267a 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -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;
}
/**
diff --git a/src/multipole.h b/src/multipole.h
index e867dfd4e2cc5c9fcd06d7d95dcf76a97689c2b3..fe17a624491ac88a8d9061cba99b47c69d48bc74 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -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]);
diff --git a/src/serial_io.c b/src/serial_io.c
index 3b3872762f237c4922bbc30fd866599547b23557..05a1839572cfcea8a81559cbdb197dbed26e605d 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -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,
- int n_threads, int dry_run) {
+ size_t* Nblackholes, int* flag_entropy, float* gpart_mass,
+ 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, int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -801,6 +801,14 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Record the mass of the DM particles */
+ const float local_gpart_mass = (*gparts)[0].mass;
+ float global_gpart_mass;
+ MPI_Allreduce(&local_gpart_mass, &global_gpart_mass, 1, MPI_FLOAT, MPI_MAX,
+ comm);
+ /* Record the mass of the DM particles */
+ *gpart_mass = global_gpart_mass;
+
/* Prepare the DM background particles */
io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
diff --git a/src/serial_io.h b/src/serial_io.h
index 6badf076b8373ce2b6d5b2cd3d91a5b2f57ed8ae..e038d578a8d8c1acfd582515263ef9b086f9edbe 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -39,11 +39,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,
- int n_threads, int dry_run);
+ size_t* Nblackholes, int* flag_entropy, float* gpart_mass,
+ 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, int n_threads, int dry_run);
void write_output_serial(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/single_io.c b/src/single_io.c
index 6f1f4a0d8bbcb1c66c35e6efe6f2a5a13a9f5104..ae42fd9a213e05ad9688e4642f0f2fb989302964 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -405,10 +405,10 @@ void read_ic_single(const char* fileName,
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 n_threads, int dry_run) {
+ size_t* Nblackholes, int* flag_entropy, float* gpart_mass,
+ 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 n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -665,6 +665,9 @@ void read_ic_single(const char* fileName,
/* Prepare the DM particles */
io_prepare_dm_gparts(&tp, *gparts, Ndm);
+ /* Record the mass of the DM particles */
+ *gpart_mass = (*gparts)[0].mass;
+
/* Prepare the DM background particles */
io_prepare_dm_background_gparts(&tp, *gparts + Ndm, Ndm_background);
diff --git a/src/single_io.h b/src/single_io.h
index 47b575272c2fbe5b170eb4370743cd6830b99ae0..6044d903efffe587d32d66e4e6785503ed6b10d2 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -35,10 +35,10 @@ void read_ic_single(const char* fileName,
struct part** parts, struct gpart** gparts,
struct spart** sparts, struct bpart** bparts, size_t* Ngas,
size_t* Ndm, size_t* Ndm_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 nr_threads, int dry_run);
+ size_t* Nblackholes, int* flag_entropy, float* gpart_mass,
+ 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 nr_threads, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/space.c b/src/space.c
index 61ecf3fd81f4a57d274ad911e0471337625f944b..cea579d09f438ff07611eb4c674decd7a3547de8 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1928,7 +1928,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
/* Check that the multipole construction went OK */
if (s->with_self_gravity)
for (int k = 0; k < s->nr_cells; k++)
- cell_check_multipole(&s->cells_top[k]);
+ cell_check_multipole(&s->cells_top[k], s->e->gravity_properties);
#endif
/* Clean up any stray sort indices in the cell buffer. */
@@ -3603,7 +3603,8 @@ void space_split_recursive(struct space *s, struct cell *c,
if (s->with_self_gravity) {
if (gcount > 0) {
- gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count);
+ gravity_P2M(c->grav.multipole, c->grav.parts, c->grav.count,
+ e->gravity_properties);
} else {