diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
index e7c2538c6238b3aaaace51e54bdc8331b5ea9b05..e96d3520a08e12228774f13d7d0740afc009669d 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml
@@ -37,7 +37,7 @@ Statistics:
InitialConditions:
file_name: lowres8.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
- stars_smoothing_length: 0.3
+ stars_smoothing_length: 0.5
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/src/runner.c b/src/runner.c
index 9ca01904bca1fdfb2a7fe1d3aa1f1228482213a0..658545806403b02b92eda2ed065d7353c6f2f40f 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -223,7 +223,9 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
#ifdef SWIFT_DEBUG_CHECKS
if ((f > 0.f && h_new > h_old) || (f < 0.f && h_new < h_old))
error(
- "Smoothing length correction not going in the right direction");
+ "Smoothing length correction not going in the right direction"
+ "sp->id=%lld sp->h=%e",
+ sp->id, sp->h);
#endif
/* Safety check: truncate to the range [ h_old/2 , 2h_old ]. */
diff --git a/src/stars/EAGLE/stars.h b/src/stars/EAGLE/stars.h
index a1817fe4947043d18312e7180ac02d3c68529f9d..324cda396d5bbfbbd9df30e62f2ddec8f2c68ac7 100644
--- a/src/stars/EAGLE/stars.h
+++ b/src/stars/EAGLE/stars.h
@@ -33,21 +33,6 @@ __attribute__((always_inline)) INLINE static float stars_compute_timestep(
return FLT_MAX;
}
-/**
- * @brief Initialises the s-particles for the first time
- *
- * This function is called only once just after the ICs have been
- * read in to do some conversions.
- *
- * @param sp The particle to act upon
- */
-__attribute__((always_inline)) INLINE static void stars_first_init_spart(
- struct spart* sp) {
-
- sp->time_bin = 0;
- sp->birth_density = -1.f;
-}
-
/**
* @brief Prepares a s-particle for its interactions
*
@@ -64,6 +49,24 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
sp->density.wcount = 0.f;
sp->density.wcount_dh = 0.f;
+ sp->rho_gas = 0.f;
+}
+
+/**
+ * @brief Initialises the s-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void stars_first_init_spart(
+ struct spart* sp) {
+
+ sp->time_bin = 0;
+ sp->birth_density = -1.f;
+
+ stars_init_spart(sp);
}
/**
@@ -132,6 +135,7 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
/* Finish the calculation by inserting the missing h-factors */
+ sp->rho_gas *= h_inv_dim;
sp->density.wcount *= h_inv_dim;
sp->density.wcount_dh *= h_inv_dim_plus_one;
}
@@ -146,14 +150,10 @@ __attribute__((always_inline)) INLINE static void stars_end_density(
__attribute__((always_inline)) INLINE static void stars_spart_has_no_neighbours(
struct spart* restrict sp, const struct cosmology* cosmo) {
- /* Some smoothing length multiples. */
- const float h = sp->h;
- const float h_inv = 1.0f / h; /* 1/h */
- const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
-
/* Re-set problematic values */
- sp->density.wcount = kernel_root * h_inv_dim;
+ sp->density.wcount = 0.f;
sp->density.wcount_dh = 0.f;
+ sp->rho_gas = 0.f;
}
/**
diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h
index 9e27f86028245a230cfd777dfc46da7b7d2f3915..5c4d94c2060ec42633a5eec472d844ca919b56eb 100644
--- a/src/stars/EAGLE/stars_iact.h
+++ b/src/stars/EAGLE/stars_iact.h
@@ -16,6 +16,9 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
const struct part *restrict pj, float a,
float H) {
+ /* Get the gas mass. */
+ const float mj = pj->mass;
+
float wi, wi_dx;
/* Get r and 1/r. */
@@ -31,6 +34,9 @@ runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj,
si->density.wcount += wi;
si->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+ /* Compute contribution to the density */
+ si->rho_gas += mj * wi;
+
#ifdef DEBUG_INTERACTIONS_STARS
/* Update ngb counters */
if (si->num_ngb_density < MAX_NUM_OF_NEIGHBOURS_STARS)
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 64fde1a779966cc072a750fab89c74c5329bc03b..8a5f0d16551afc01a91a211f20198ed33ea262c6 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -62,7 +62,7 @@ INLINE static void stars_write_particles(const struct spart *sparts,
int *num_fields) {
/* Say how much we want to write */
- *num_fields = 8;
+ *num_fields = 9;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -77,10 +77,12 @@ INLINE static void stars_write_particles(const struct spart *sparts,
sparts, h);
list[5] = io_make_output_field("BirthDensity", FLOAT, 1, UNIT_CONV_DENSITY,
sparts, birth_density);
- list[6] = io_make_output_field("Initial_Masses", FLOAT, 1, UNIT_CONV_MASS,
+ list[6] = io_make_output_field("InitialMasses", FLOAT, 1, UNIT_CONV_MASS,
sparts, mass_init);
- list[7] = io_make_output_field("Birth_time", FLOAT, 1, UNIT_CONV_TIME, sparts,
+ list[7] = io_make_output_field("Birthtime", FLOAT, 1, UNIT_CONV_TIME, sparts,
birth_time);
+ list[8] = io_make_output_field("GasDensity", FLOAT, 1, UNIT_CONV_DENSITY,
+ sparts, rho_gas);
}
/**
diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h
index 02ba397660302a06748b3db9c4b3371332c0478a..718003c4080913b65ad64f002f5e23a4f031c2dc 100644
--- a/src/stars/EAGLE/stars_part.h
+++ b/src/stars/EAGLE/stars_part.h
@@ -58,9 +58,12 @@ struct spart {
/*! Initial star mass */
float mass_init;
- /* Particle cutoff radius. */
+ /*! Particle smoothing length. */
float h;
+ /*! Density of the gas surrounding the star. */
+ float rho_gas;
+
/*! Particle time bin */
timebin_t time_bin;