Added the option to generate gas particles directly from the ICs.

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);

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.

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.;

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);

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 {

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
  *

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);