Can correctly read/write star particles alongside other combinations of gas/DM

examples/MultiTypes/makeIC.py

@@ -36,6 +36,9 @@ eta = 1.2349             # 48 ngbs with cubic spline kernel
 rhoDM = 1.
 Ldm = int(sys.argv[2])  # Number of particles along one axis
 
+massStars = 0.1
+Lstars = int(sys.argv[3])  # Number of particles along one axis
+
 fileName = "multiTypes.hdf5"
 
 #---------------------------------------------------
@@ -46,6 +49,10 @@ internalEnergy = P / ((gamma - 1.)*rhoGas)
 numDM = Ldm**3
 massDM = boxSize**3 * rhoDM / numDM
 
+numStars = Lstars**3
+massStars = massDM * massStars
+
+
 #--------------------------------------------------
 
 #File
@@ -54,9 +61,9 @@ file = h5py.File(fileName, 'w')
 # Header
 grp = file.create_group("/Header")
 grp.attrs["BoxSize"] = boxSize
-grp.attrs["NumPart_Total"] =  [numGas, numDM, 0, 0, 0, 0]
+grp.attrs["NumPart_Total"] =  [numGas, numDM, 0, 0, numStars, 0]
 grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
-grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, numStars, 0]
 grp.attrs["Time"] = 0.0
 grp.attrs["NumFilesPerSnapshot"] = 1
 grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
@@ -142,4 +149,33 @@ coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
 ds = grp.create_dataset('Coordinates', (numDM, 3), 'd')
 ds[()] = coords
 
+
+
+# Star Particle group
+grp = file.create_group("/PartType4")
+
+v  = zeros((numStars, 3))
+ds = grp.create_dataset('Velocities', (numStars, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numStars, 1), massStars)
+ds = grp.create_dataset('Masses', (numStars,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+ids = linspace(0, numStars, numStars, endpoint=False).reshape((numStars,1))
+ds = grp.create_dataset('ParticleIDs', (numStars, 1), 'L')
+ds[()] = ids + Lgas**3 + 1
+x      = ids % Ldm;
+y      = ((ids - x) / Ldm) % Ldm;
+z      = (ids - x - Ldm * y) / Ldm**2;
+coords = zeros((numStars, 3))
+coords[:,0] = z[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+coords[:,1] = y[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+coords[:,2] = x[:,0] * boxSize / Ldm + boxSize / (2*Ldm)
+ds = grp.create_dataset('Coordinates', (numStars, 3), 'd')
+ds[()] = coords
+
+
 file.close()

examples/main.c

@@ -393,7 +393,9 @@ int main(int argc, char *argv[]) {
 #endif
 #else
   read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
-                 &Nspart, &periodic, &flag_entropy_ICs, dry_run);
+                 &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+                 (with_external_gravity || with_self_gravity), with_stars,
+                 dry_run);
 #endif
   if (myrank == 0) {
     clocks_gettime(&toc);
@@ -414,7 +416,7 @@ int main(int argc, char *argv[]) {
   if (!with_stars) {
     free(sparts);
     sparts = NULL;
-    for (size_t k = 0; k < Nspart; ++k)
+    for (size_t k = 0; k < Ngpart; ++k)
       if (gparts[k].type == swift_type_star) error("Linking problem");
     Nspart = 0;
   }

src/single_io.c

@@ -337,6 +337,7 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
                     struct spart** sparts, size_t* Ngas, size_t* Ngparts,
                     size_t* Nstars, int* periodic, int* flag_entropy,
+                    int with_hydro, int with_gravity, int with_stars,
                     int dry_run) {
 
   hid_t h_file = 0, h_grp = 0;
@@ -443,26 +444,32 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
         units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
 
   /* Allocate memory to store SPH particles */
+  if (with_hydro) {
     *Ngas = N[GAS];
     if (posix_memalign((void*)parts, part_align, *Ngas * sizeof(struct part)) !=
         0)
       error("Error while allocating memory for SPH particles");
     bzero(*parts, *Ngas * sizeof(struct part));
+  }
 
   /* Allocate memory to store star particles */
+  if (with_stars) {
     *Nstars = N[STAR];
     if (posix_memalign((void*)sparts, spart_align,
                        *Nstars * sizeof(struct spart)) != 0)
       error("Error while allocating memory for star particles");
     bzero(*sparts, *Nstars * sizeof(struct spart));
+  }
 
-  /* Allocate memory to store all particles */
+  /* Allocate memory to store all gravity particles */
+  if (with_gravity) {
     Ndm = N[DM];
-  *Ngparts = N[GAS] + N[DM] + N[STAR];
+    *Ngparts = (with_hydro ? N[GAS] : 0) + N[DM] + (with_stars ? N[STAR] : 0);
     if (posix_memalign((void*)gparts, gpart_align,
                        *Ngparts * sizeof(struct gpart)) != 0)
       error("Error while allocating memory for gravity particles");
     bzero(*gparts, *Ngparts * sizeof(struct gpart));
+  }
 
   /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
    * (1024.*1024.)); */
@@ -493,18 +500,24 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
     switch (ptype) {
 
       case GAS:
+        if (with_hydro) {
           Nparticles = *Ngas;
           hydro_read_particles(*parts, list, &num_fields);
+        }
         break;
 
       case DM:
+        if (with_gravity) {
           Nparticles = Ndm;
           darkmatter_read_particles(*gparts, list, &num_fields);
+        }
         break;
 
       case STAR:
+        if (with_stars) {
           Nparticles = *Nstars;
           star_read_particles(*sparts, list, &num_fields);
+        }
         break;
 
       default:
@@ -521,13 +534,15 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
   }
 
   /* Prepare the DM particles */
-  if (!dry_run) prepare_dm_gparts(*gparts, Ndm);
+  if (!dry_run && with_gravity) prepare_dm_gparts(*gparts, Ndm);
 
   /* Duplicate the hydro particles into gparts */
-  if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+  if (!dry_run && with_gravity && with_hydro)
+    duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
 
   /* Duplicate the star particles into gparts */
-  if (!dry_run) duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
+  if (!dry_run && with_gravity && with_stars)
+    duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
 
   /* message("Done Reading particles..."); */
 
@@ -765,7 +780,10 @@ void write_output_single(struct engine* e, const char* baseName,
                  internal_units, snapshot_units);
 
     /* Free temporary array */
+    if (dmparts) {
       free(dmparts);
+      dmparts = NULL;
+    }
 
     /* Close particle group */
     H5Gclose(h_grp);

src/single_io.h

@@ -33,6 +33,7 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
                     double dim[3], struct part** parts, struct gpart** gparts,
                     struct spart** sparts, size_t* Ngas, size_t* Ndm,
                     size_t* Nstars, int* periodic, int* flag_entropy,
+                    int with_hydro, int with_gravity, int with_stars,
                     int dry_run);
 
 void write_output_single(struct engine* e, const char* baseName,