diff --git a/examples/main.c b/examples/main.c
index 90528e94fdb5c07417a68d212cdb7d9b0b2ddcb6..9c9370a2b1e14b09563aca7f3dc017d7041a336c 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -624,6 +624,12 @@ int main(int argc, char *argv[]) {
params, "InitialConditions:cleanup_smoothing_lengths", 0);
const int cleanup_h = parser_get_opt_param_int(
params, "InitialConditions:cleanup_h_factors", 0);
+ const int generate_gas_in_ics = parser_get_opt_param_int(
+ params, "InitialConditions:generate_gas_in_ics", 0);
+ 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.");
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
if (myrank == 0 && cleanup_h)
message("Cleaning up h-factors (h=%f)", cosmo.h);
@@ -694,8 +700,9 @@ int main(int argc, char *argv[]) {
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
- space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
- periodic, replicate, with_self_gravity, talking, dry_run);
+ space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
+ Nspart, periodic, replicate, generate_gas_in_ics,
+ with_self_gravity, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 07c02b76bfd36cd9021e3420d3fce9a1383b4fb4..ff70566644f7334145dd35adbe264f4de8cc56d2 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -81,6 +81,7 @@ Statistics:
# Parameters related to the initial conditions
InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read
+ generate_gas_in_ics: 0 # (Optional) Generate gas particles from the DM-only ICs (e.g. from panphasia).
cleanup_h_factors: 0 # (Optional) Clean up the h-factors used in the ICs (e.g. in Gadget files).
cleanup_smoothing_lengths: 0 # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
smoothing_length_scaling: 1. # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
diff --git a/src/cosmology.c b/src/cosmology.c
index 15a2d09f909f5de2675ec8de12ab400c92d18beb..23fcddcf61e8082568b24130dc807565a4ab7004 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -130,9 +130,11 @@ double cosmology_get_time_since_big_bang(const struct cosmology *c, double a) {
* @brief Update the cosmological parameters to the current simulation time.
*
* @param c The #cosmology struct.
+ * @param phys_const The physical constants in the internal units.
* @param ti_current The current (integer) time.
*/
-void cosmology_update(struct cosmology *c, integertime_t ti_current) {
+void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
+ integertime_t ti_current) {
/* Get scale factor and powers of it */
const double a = c->a_begin * exp(ti_current * c->time_base);
@@ -173,6 +175,10 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) {
/* Expansion rate */
c->a_dot = c->H * c->a;
+ /* Critical density */
+ c->critical_density =
+ 3. * c->H * c->H / (8. * M_PI * phys_const->const_newton_G);
+
/* Time-step conversion factor */
c->time_step_factor = c->H;
@@ -431,7 +437,7 @@ void cosmology_init(const struct swift_params *params,
cosmology_init_tables(c);
/* Set remaining variables to alid values */
- cosmology_update(c, 0);
+ cosmology_update(c, phys_const, 0);
/* Update the times */
c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
@@ -474,6 +480,8 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.;
+ c->critical_density = 0.;
+
c->time_step_factor = 1.;
c->a_dot = 0.;
diff --git a/src/cosmology.h b/src/cosmology.h
index b541e20fa1dd4f55422f9ff1d0952783d8de0ad7..109b80a57d8dbc4eb942dd4ecbbc0db84198100b 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -70,6 +70,9 @@ struct cosmology {
/*! Hubble constant at the current redshift (in internal units) */
double H;
+ /*! The critical density at the current redshift (in internal units) */
+ double critical_density;
+
/*! Conversion factor from internal time-step size to cosmological step */
double time_step_factor;
@@ -160,7 +163,8 @@ struct cosmology {
double universe_age_at_present_day;
};
-void cosmology_update(struct cosmology *c, integertime_t ti_current);
+void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
+ integertime_t ti_current);
double cosmology_get_drift_factor(const struct cosmology *cosmo,
integertime_t ti_start, integertime_t ti_end);
diff --git a/src/engine.c b/src/engine.c
index 9abb75936052c1f5fad6c4efe6170467a3111a49..962aa0e2108994371bc9dd0c16b2b07f6e6be122 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4395,7 +4395,7 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_cosmology) {
e->time_old = e->time;
- cosmology_update(e->cosmology, e->ti_current);
+ cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
e->time = e->cosmology->time;
e->time_step = e->time - e->time_old;
} else {
diff --git a/src/space.c b/src/space.c
index fa080f784443fb801d963e08ea7f1ec42e7b7caf..51107706684be5d3cbdd959597214d541240d3d5 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2854,7 +2854,8 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param Ngpart The number of Gravity particles in the space.
* @param Nspart The number of star particles in the space.
* @param periodic flag whether the domain is periodic or not.
- * @param replicate How many replications along each direction do we want ?
+ * @param replicate How many replications along each direction do we want?
+ * @param generate_gas_in_ics Are we generating gas particles from the gparts?
* @param gravity flag whether we are doing gravity or not.
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
@@ -2865,9 +2866,10 @@ void space_convert_quantities(struct space *s, int verbose) {
* recursively.
*/
void space_init(struct space *s, const struct swift_params *params,
- double dim[3], struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t Npart, size_t Ngpart,
- size_t Nspart, int periodic, int replicate, int gravity,
+ const struct cosmology *cosmo, double dim[3],
+ struct part *parts, struct gpart *gparts, struct spart *sparts,
+ size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
+ int replicate, int generate_gas_in_ics, int self_gravity,
int verbose, int dry_run) {
/* Clean-up everything */
@@ -2878,7 +2880,7 @@ void space_init(struct space *s, const struct swift_params *params,
s->dim[1] = dim[1];
s->dim[2] = dim[2];
s->periodic = periodic;
- s->gravity = gravity;
+ s->gravity = self_gravity;
s->nr_parts = Npart;
s->size_parts = Npart;
s->parts = parts;
@@ -2890,6 +2892,9 @@ void space_init(struct space *s, const struct swift_params *params,
s->sparts = sparts;
s->nr_queues = 1; /* Temporary value until engine construction */
+ /* Are we generating gas from the DM-only ICs? */
+ if (generate_gas_in_ics) space_generate_gas(s, cosmo, verbose);
+
/* Are we replicating the space ? */
if (replicate < 1)
error("Value of 'InitialConditions:replicate' (%d) is too small",
@@ -3166,6 +3171,91 @@ void space_replicate(struct space *s, int replicate, int verbose) {
#endif
}
+void space_generate_gas(struct space *s, const struct cosmology *cosmo,
+ int verbose) {
+
+ if (verbose) message("Generating gas particles from gparts");
+
+ /* Store the current values */
+ const size_t nr_parts = s->nr_parts;
+ const size_t nr_gparts = s->nr_gparts;
+
+ if (nr_parts != 0)
+ error("Generating gas particles from DM but gas already exists!");
+
+ s->size_parts = s->nr_parts = nr_gparts;
+ s->size_gparts = s->nr_gparts = 2 * nr_gparts;
+
+ /* Allocate space for new particles */
+ struct part *parts = NULL;
+ struct gpart *gparts = NULL;
+
+ if (posix_memalign((void **)&parts, part_align,
+ s->nr_parts * sizeof(struct part)) != 0)
+ error("Failed to allocate new part array.");
+
+ if (posix_memalign((void **)&gparts, gpart_align,
+ s->nr_gparts * sizeof(struct gpart)) != 0)
+ error("Failed to allocate new gpart array.");
+
+ /* Start by copying the gparts */
+ memcpy(gparts + nr_gparts, gparts, nr_gparts * sizeof(struct gpart));
+
+ /* And zero the parts */
+ bzero(parts, s->nr_parts * sizeof(struct part));
+
+ /* Compute some constants */
+ const double mass_ratio = cosmo->Omega_b / cosmo->Omega_m;
+ const double bg_density = cosmo->Omega_m * cosmo->critical_density;
+ const double bg_density_inv = 1. / bg_density;
+
+ /* Update the particle properties */
+ for (size_t i = 0; i < nr_gparts; ++i) {
+
+ struct part *p = &parts[i];
+ struct gpart *gp_dm = &gparts[i];
+ struct gpart *gp_gas = &gparts[nr_gparts + i];
+
+ /* Set the IDs */
+ p->id = gp_gas->id_or_neg_offset * 2 + 1;
+ gp_dm->id_or_neg_offset *= 2;
+
+ /* Set the links correctly */
+ p->gpart = gp_gas;
+ gp_gas->id_or_neg_offset = -i;
+ gp_gas->type = swift_type_gas;
+
+ /* Compute positions shift */
+ const double d = cbrt(gp_dm->mass * bg_density_inv);
+ const double shift_dm = d * mass_ratio;
+ const double shift_gas = d * (1. - mass_ratio);
+
+ /* Set the masses */
+ gp_dm->mass *= (1. - mass_ratio);
+ gp_gas->mass *= mass_ratio;
+ p->mass = gp_gas->mass;
+
+ /* Set the new positions */
+ gp_dm->x[0] += shift_dm;
+ gp_dm->x[1] += shift_dm;
+ gp_dm->x[2] += shift_dm;
+ gp_gas->x[0] += shift_gas;
+ gp_gas->x[1] += shift_gas;
+ gp_gas->x[2] += shift_gas;
+ p->x[0] = gp_gas->x[0];
+ p->x[1] = gp_gas->x[1];
+ p->x[2] = gp_gas->x[2];
+
+ /* Also copy the velocities */
+ p->v[0] = gp_gas->v_full[0];
+ p->v[1] = gp_gas->v_full[1];
+ p->v[2] = gp_gas->v_full[2];
+
+ /* Set the smoothing length to the mean inter-particle separation */
+ p->h = d;
+ }
+}
+
/**
* @brief Cleans-up all the cell links in the space
*
diff --git a/src/space.h b/src/space.h
index 11cbaabdc9fc3ae17024c042ae868d464954d501..76ff2ea9908195d618663eb1874eb48bec463bb7 100644
--- a/src/space.h
+++ b/src/space.h
@@ -38,6 +38,7 @@
/* Avoid cyclic inclusions */
struct cell;
+struct cosmology;
/* Some constants. */
#define space_cellallocchunk 1000
@@ -182,9 +183,10 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
void space_init(struct space *s, const struct swift_params *params,
- double dim[3], struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t Npart, size_t Ngpart,
- size_t Nspart, int periodic, int replicate, int gravity,
+ const struct cosmology *cosmo, double dim[3],
+ struct part *parts, struct gpart *gparts, struct spart *sparts,
+ size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
+ int replicate, int generate_gas_in_ics, int gravity,
int verbose, int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
@@ -239,6 +241,8 @@ void space_check_top_multipoles_drift_point(struct space *s,
integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
+void space_generate_gas(struct space *s, const struct cosmology *cosmo,
+ int verbose);
void space_reset_task_counters(struct space *s);
void space_clean(struct space *s);
void space_free_cells(struct space *s);