diff --git a/examples/main.c b/examples/main.c
index 38847543bd25d6ce877b713b7e6e9ff614eb16bb..2542be95e0b02465fa5b1bdac9b2c2757633a787 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -135,6 +135,7 @@ int main(int argc, char *argv[]) {
struct gpart *gparts = NULL;
struct gravity_props gravity_properties;
struct hydro_props hydro_properties;
+ struct stars_props stars_properties;
struct part *parts = NULL;
struct phys_const prog_const;
struct sourceterms sourceterms;
@@ -673,6 +674,12 @@ int main(int argc, char *argv[]) {
else
bzero(&eos, sizeof(struct eos_parameters));
+ /* Initialise the stars properties */
+ if (with_stars)
+ stars_props_init(&stars_properties, &prog_const, &us, params, &hydro_properties);
+ else
+ bzero(&stars_properties, sizeof(struct stars_props));
+
/* Initialise the gravity properties */
if (with_self_gravity)
gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology);
@@ -886,7 +893,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
engine_init(&e, &s, params, N_total[0], N_total[1], N_total[2],
engine_policies, talking, &reparttype, &us, &prog_const, &cosmo,
- &hydro_properties, &gravity_properties, &mesh, &potential,
+ &hydro_properties, &gravity_properties, &stars_properties, &mesh, &potential,
&cooling_func, &chemistry, &sourceterms);
engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff,
talking, restart_file);
diff --git a/src/engine.c b/src/engine.c
index b79e1296d51478ca5569dce5fbab986d306d1ba4..adb853ac91d89c48855ba2b46fc058a5a7805d40 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -5776,6 +5776,7 @@ void engine_unpin(void) {
* @param cosmo The #cosmology used for this run.
* @param hydro The #hydro_props used for this run.
* @param gravity The #gravity_props used for this run.
+ * @param stars The #stars_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.
* @param cooling_func The properties of the cooling function.
@@ -5788,7 +5789,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
const struct unit_system *internal_units,
const struct phys_const *physical_constants,
struct cosmology *cosmo, const struct hydro_props *hydro,
- struct gravity_props *gravity, struct pm_mesh *mesh,
+ struct gravity_props *gravity, const struct stars_props *stars,
+ struct pm_mesh *mesh,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_global_data *chemistry,
@@ -5854,6 +5856,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->cosmology = cosmo;
e->hydro_properties = hydro;
e->gravity_properties = gravity;
+ e->stars_properties = stars;
e->mesh = mesh;
e->external_potential = potential;
e->cooling_func = cooling_func;
@@ -6187,6 +6190,9 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
if (e->policy & engine_policy_self_gravity)
if (e->nodeID == 0) gravity_props_print(e->gravity_properties);
+ if (e->policy & engine_policy_stars)
+ if (e->nodeID == 0) stars_props_print(e->stars_properties);
+
/* Check we have sensible time bounds */
if (e->time_begin >= e->time_end)
error(
@@ -6915,6 +6921,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
phys_const_struct_dump(e->physical_constants, stream);
hydro_props_struct_dump(e->hydro_properties, stream);
gravity_props_struct_dump(e->gravity_properties, stream);
+ stars_props_struct_dump(e->stars_properties, stream);
pm_mesh_struct_dump(e->mesh, stream);
potential_struct_dump(e->external_potential, stream);
cooling_struct_dump(e->cooling_func, stream);
@@ -6990,6 +6997,11 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
gravity_props_struct_restore(gravity_properties, stream);
e->gravity_properties = gravity_properties;
+ struct stars_props *stars_properties =
+ (struct stars_props *)malloc(sizeof(struct stars_props));
+ stars_props_struct_restore(stars_properties, stream);
+ e->stars_properties = stars_properties;
+
struct pm_mesh *mesh = (struct pm_mesh *)malloc(sizeof(struct pm_mesh));
pm_mesh_struct_restore(mesh, stream);
e->mesh = mesh;
diff --git a/src/engine.h b/src/engine.h
index 65cb7f22a715d01c745268948c796db3a1f735ab..5bb2780211f953edefa192e90e4ae1076409c4a5 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -323,6 +323,9 @@ struct engine {
/* Properties of the hydro scheme */
const struct hydro_props *hydro_properties;
+ /* Properties of the star model */
+ const struct stars_props *stars_properties;
+
/* Properties of the self-gravity scheme */
struct gravity_props *gravity_properties;
diff --git a/src/runner.c b/src/runner.c
index b1b783847fd141f39c47ca6362116ffd1549a6cf..cbd4e6d7535037356997b286d256b740697388f9 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -149,11 +149,12 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct space *s = e->s;
const struct cosmology *cosmo = e->cosmology;
- const float star_h_max = e->star_properties->h_max;
- const float eps = e->star_properties->h_tolerance;
+ const struct stars_props *stars_properties = e->stars_properties;
+ const float star_h_max = e->stars_properties->h_max;
+ const float eps = e->stars_properties->h_tolerance;
const float star_eta_dim =
- pow_dimension(e->star_properties->eta_neighbours);
- const int max_smoothing_iter = e->star_properties->max_smoothing_iterations;
+ pow_dimension(e->stars_properties->eta_neighbours);
+ const int max_smoothing_iter = e->stars_properties->max_smoothing_iterations;
int redo = 0, count = 0;
TIMER_TIC;
@@ -212,7 +213,7 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
} else {
/* Finish the density calculation */
- star_end_density(p, cosmo);
+ star_end_density(sp, cosmo);
/* Compute one step of the Newton-Raphson scheme */
const float n_sum = sp->wcount * h_old_dim;
@@ -250,7 +251,7 @@ void runner_do_star_ghost(struct runner *r, struct cell *c, int timer) {
redo += 1;
/* Re-initialise everything */
- star_init_spart(sp, hs);
+ star_init_spart(sp, stars_properties);
/* Off we go ! */
continue;
diff --git a/src/runner_doiact_star.h b/src/runner_doiact_star.h
index f59069935e778b4c72449499eaf2e0bdbc89e523..555b665c81b88380c9c1a1a96d292c521a203ed1 100644
--- a/src/runner_doiact_star.h
+++ b/src/runner_doiact_star.h
@@ -155,7 +155,7 @@ void runner_dosubpair_star_density(struct runner *r, struct cell *restrict ci,
#endif
if (r2 < hig2)
- runner_iact_nonsym_star_density(r2, dx, hj, hi, si, pj, a, H);
+ runner_iact_nonsym_star_density(r2, dx, hj, hi, si, pj, a, H);
} /* loop over the parts in cj. */
} /* loop over the parts in ci. */
@@ -173,3 +173,81 @@ void runner_dopair_star_density(struct runner *r, struct cell *restrict ci,
if (timer) TIMER_TOC(timer_dopair_star_density);
}
+
+
+/**
+ * @brief Compute the interactions between a cell pair, but only for the
+ * given indices in ci.
+ *
+ * Version using a brute-force algorithm.
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param sparts_i The #part to interact with @c cj.
+ * @param ind The list of indices of particles in @c ci to interact with.
+ * @param scount The number of particles in @c ind.
+ * @param cj The second #cell.
+ * @param shift The shift vector to apply to the particles in ci.
+ */
+void runner_dopair_subset_star_density(struct runner *r, struct cell *restrict ci,
+ struct spart *restrict sparts_i, int *restrict ind,
+ int scount, struct cell *restrict cj,
+ const double *shift) {
+
+ const struct engine *e = r->e;
+ const struct cosmology *cosmo = e->cosmology;
+
+ TIMER_TIC;
+
+ const int count_j = cj->count;
+ struct part *restrict parts_j = cj->parts;
+
+ /* Cosmological terms */
+ const float a = cosmo->a;
+ const float H = cosmo->H;
+
+ /* Loop over the parts_i. */
+ for (int pid = 0; pid < scount; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct spart *restrict spi = &sparts_i[ind[pid]];
+ double spix[3];
+ for (int k = 0; k < 3; k++) spix[k] = spi->x[k] - shift[k];
+ const float hi = spi->h;
+ const float hig2 = hi * hi * kernel_gamma2;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!spart_is_active(spi, e))
+ error("Trying to correct smoothing length of inactive particle !");
+#endif
+
+ /* Loop over the parts in cj. */
+ for (int pjd = 0; pjd < count_j; pjd++) {
+
+ /* Get a pointer to the jth particle. */
+ struct part *restrict pj = &parts_j[pjd];
+
+ /* Compute the pairwise distance. */
+ float r2 = 0.0f;
+ float dx[3];
+ for (int k = 0; k < 3; k++) {
+ dx[k] = spix[k] - pj->x[k];
+ r2 += dx[k] * dx[k];
+ }
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (spi->ti_drift != e->ti_current)
+ error("Particle pi not drifted to current time");
+ if (pj->ti_drift != e->ti_current)
+ error("Particle pj not drifted to current time");
+#endif
+ /* Hit or miss? */
+ if (r2 < hig2) {
+ runner_iact_nonsym_star_density(r2, dx, hi, pj->h, spi, pj, a, H);
+ }
+ } /* loop over the parts in cj. */
+ } /* loop over the parts in ci. */
+
+ TIMER_TOC(timer_dopair_subset_naive);
+}
diff --git a/src/stars/Default/star.h b/src/stars/Default/star.h
index 41ee4f7887a92838f39803dc348c653e406273ad..325cd26b3883d59b90341abd4ff41204ac789b3a 100644
--- a/src/stars/Default/star.h
+++ b/src/stars/Default/star.h
@@ -53,7 +53,7 @@ __attribute__((always_inline)) INLINE static void star_first_init_spart(
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void star_init_spart(
- struct spart* sp) {
+ struct spart* sp, const struct stars_props *stars_properties) {
#ifdef DEBUG_INTERACTIONS_SPH
for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS; ++i) sp->ids_ngbs_density[i] = -1;
@@ -98,7 +98,7 @@ __attribute__((always_inline)) INLINE static void star_kick_extra(
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void star_end_density(
- struct spart* sp) {
+ struct spart* sp, const struct cosmology* cosmo) {
/* Some smoothing length multiples. */
const float h = sp->h;
@@ -143,7 +143,7 @@ __attribute__((always_inline)) INLINE static void star_spart_has_no_neighbours(
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void star_prepare_force(
- struct part *restrict sp, const struct cosmology *cosmo) {}
+ struct spart *restrict sp, const struct cosmology *cosmo) {}
/**
diff --git a/src/stars/Default/star_io.h b/src/stars/Default/star_io.h
index 7ad29f0a935c002b1337c2a75d6f987c05c9bb43..3c424dc0ae703896605ccb8f23b7b9684f0838b9 100644
--- a/src/stars/Default/star_io.h
+++ b/src/stars/Default/star_io.h
@@ -19,6 +19,7 @@
#ifndef SWIFT_DEFAULT_STAR_IO_H
#define SWIFT_DEFAULT_STAR_IO_H
+#include "star_part.h"
#include "io_properties.h"
/**
@@ -70,4 +71,117 @@ INLINE static void star_write_particles(const struct spart* sparts,
sparts, id);
}
+/**
+ * @brief Initialize the global properties of the stars scheme.
+ *
+ * @param p The #stars_props.
+ * @param phys_const The physical constants in the internal unit system.
+ * @param us The internal unit system.
+ * @param params The parsed parameters.
+ */
+INLINE static void stars_props_init(struct stars_props *sp,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ struct swift_params *params,
+ const struct hydro_props *p) {
+
+ /* Kernel properties */
+ sp->eta_neighbours = parser_get_param_float(params, "Stars:resolution_eta");
+
+ /* Tolerance for the smoothing length Newton-Raphson scheme */
+ sp->h_tolerance = parser_get_opt_param_float(params, "Stars:h_tolerance",
+ p->h_tolerance);
+
+ /* Get derived properties */
+ sp->target_neighbours = pow_dimension(sp->eta_neighbours) * kernel_norm;
+ const float delta_eta = sp->eta_neighbours * (1.f + sp->h_tolerance);
+ sp->delta_neighbours =
+ (pow_dimension(delta_eta) - pow_dimension(sp->eta_neighbours)) *
+ kernel_norm;
+
+ /* Maximal smoothing length */
+ sp->h_max = parser_get_opt_param_float(params, "Stars:h_max",
+ sp->h_max);
+
+ /* Number of iterations to converge h */
+ sp->max_smoothing_iterations = parser_get_opt_param_int(
+ params, "Stars:max_ghost_iterations", p->max_smoothing_iterations);
+
+ /* Time integration properties */
+ const float max_volume_change = parser_get_opt_param_float(
+ params, "Stars:max_volume_change", -1);
+ if (max_volume_change == -1)
+ sp->log_max_h_change = p->log_max_h_change;
+ else
+ sp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+
+}
+
+/**
+ * @brief Print the global properties of the stars scheme.
+ *
+ * @param p The #stars_props.
+ */
+void stars_props_print(const struct stars_props *sp) {
+
+ /* Now stars */
+ message("Stars kernel: %s with eta=%f (%.2f neighbours).", kernel_name,
+ sp->eta_neighbours, sp->target_neighbours);
+
+ message("Stars relative tolerance in h: %.5f (+/- %.4f neighbours).",
+ sp->h_tolerance, sp->delta_neighbours);
+
+ message(
+ "Stars integration: Max change of volume: %.2f "
+ "(max|dlog(h)/dt|=%f).",
+ pow_dimension(expf(sp->log_max_h_change)), sp->log_max_h_change);
+
+ message("Maximal smoothing length allowed: %.4f", sp->h_max);
+
+ message("Maximal iterations in ghost task set to %d",
+ sp->max_smoothing_iterations);
+
+}
+
+#if defined(HAVE_HDF5)
+void stars_props_print_snapshot(hid_t h_grpstars, const struct stars_props *sp) {
+
+ io_write_attribute_s(h_grpstars, "Kernel function", kernel_name);
+ io_write_attribute_f(h_grpstars, "Kernel target N_ngb", sp->target_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel delta N_ngb", sp->delta_neighbours);
+ io_write_attribute_f(h_grpstars, "Kernel eta", sp->eta_neighbours);
+ io_write_attribute_f(h_grpstars, "Smoothing length tolerance", sp->h_tolerance);
+ io_write_attribute_f(h_grpstars, "Maximal smoothing length", sp->h_max);
+ io_write_attribute_f(h_grpstars, "Volume log(max(delta h))",
+ sp->log_max_h_change);
+ io_write_attribute_f(h_grpstars, "Volume max change time-step",
+ pow_dimension(expf(sp->log_max_h_change)));
+ io_write_attribute_i(h_grpstars, "Max ghost iterations",
+ sp->max_smoothing_iterations);
+}
+#endif
+
+/**
+ * @brief Write a #stars_props struct to the given FILE as a stream of bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void stars_props_struct_dump(const struct stars_props *p, FILE *stream) {
+ restart_write_blocks((void *)p, sizeof(struct stars_props), 1, stream,
+ "starsprops", "stars props");
+}
+
+/**
+ * @brief Restore a stars_props struct from the given FILE as a stream of
+ * bytes.
+ *
+ * @param p the struct
+ * @param stream the file stream
+ */
+void stars_props_struct_restore(const struct stars_props *p, FILE *stream) {
+ restart_read_blocks((void *)p, sizeof(struct stars_props), 1, stream, NULL,
+ "stars props");
+}
+
#endif /* SWIFT_DEFAULT_STAR_IO_H */
diff --git a/src/stars/Default/star_part.h b/src/stars/Default/star_part.h
index e8f5ce042559db54d9c84045479353a9c5ea43d1..dd5febfb746547b7f7760cb35307a1ca5e313fce 100644
--- a/src/stars/Default/star_part.h
+++ b/src/stars/Default/star_part.h
@@ -76,4 +76,32 @@ struct spart {
} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Contains all the constants and parameters of the stars scheme
+ */
+struct stars_props {
+
+ /*! Resolution parameter */
+ float eta_neighbours;
+
+ /*! Target weightd number of neighbours (for info only)*/
+ float target_neighbours;
+
+ /*! Smoothing length tolerance */
+ float h_tolerance;
+
+ /*! Tolerance on neighbour number (for info only)*/
+ float delta_neighbours;
+
+ /*! Maximal smoothing length */
+ float h_max;
+
+ /*! Maximal number of iterations to converge h */
+ int max_smoothing_iterations;
+
+ /*! Maximal change of h over one time-step */
+ float log_max_h_change;
+};
+
#endif /* SWIFT_DEFAULT_STAR_PART_H */