diff --git a/examples/main.c b/examples/main.c
index 63753bdfceea3bb96f8e0a8b5bf36b13750335a5..f1fe32041eb5d8c3a5123b21d113befd4ab4a9e8 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -384,9 +384,10 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
- read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
- &periodic, &flag_entropy_ICs, myrank, nr_nodes,
- MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
+ read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
+ &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
+ (with_external_gravity || with_self_gravity), with_stars,
+ myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas, &Ngpart,
&Nspart, &periodic, &flag_entropy_ICs, with_hydro,
diff --git a/src/parallel_io.c b/src/parallel_io.c
index ce46d8b2dde76e28be468ac89ba0e0bf74edf786..1dbe196c151345a09ccc736958e11c6c73ea33eb 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -46,6 +46,7 @@
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
+#include "stars_io.h"
#include "units.h"
/**
@@ -373,9 +374,12 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
*/
void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
- size_t* Ngas, size_t* Ngparts, int* periodic,
- int* flag_entropy, int mpi_rank, int mpi_size,
- MPI_Comm comm, MPI_Info info, int dry_run) {
+ 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 mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
+ int dry_run) {
+
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
@@ -385,6 +389,7 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
long long N_total[NUM_PARTICLE_TYPES] = {0};
long long offset[NUM_PARTICLE_TYPES] = {0};
int dimension = 3; /* Assume 3D if nothing is specified */
+ size_t Ndm;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -492,29 +497,38 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
/* Allocate memory to store SPH particles */
- *Ngas = N[0];
- if (posix_memalign((void*)parts, part_align, (*Ngas) * sizeof(struct part)) !=
- 0)
- error("Error while allocating memory for particles");
- bzero(*parts, *Ngas * sizeof(struct part));
-
- /* Allocate memory to store all particles */
- const size_t Ndm = N[1];
- *Ngparts = N[1] + N[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) /
+ if (with_hydro) {
+ *Ngas = N[0];
+ if (posix_memalign((void*)parts, part_align, (*Ngas) * sizeof(struct part)) !=
+ 0)
+ error("Error while allocating memory for 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 gravity particles */
+ if (with_gravity) {
+ Ndm = N[1];
+ *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.)); */
/* message("BoxSize = %lf", dim[0]); */
- /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm,
- * *Ngparts); */
+ /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */
/* Loop over all particle types */
for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
@@ -539,15 +553,26 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
switch (ptype) {
case GAS:
- Nparticles = *Ngas;
- hydro_read_particles(*parts, list, &num_fields);
+ if (with_hydro) {
+ Nparticles = *Ngas;
+ hydro_read_particles(*parts, list, &num_fields);
+ }
break;
case DM:
- Nparticles = Ndm;
- darkmatter_read_particles(*gparts, list, &num_fields);
+ 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:
message("Particle Type %d not yet supported. Particles ignored", ptype);
}
@@ -563,10 +588,15 @@ void read_ic_parallel(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);
- /* Now duplicate the hydro particle into gparts */
- if (!dry_run) duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+ /* Duplicate the hydro particles into gparts */
+ if (!dry_run && with_gravity && with_hydro)
+ duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm);
+
+ /* Duplicate the star particles into gparts */
+ if (!dry_run && with_gravity && with_stars)
+ duplicate_star_gparts(*sparts, *gparts, *Nstars, Ndm + *Ngas);
/* message("Done Reading particles..."); */
@@ -609,17 +639,19 @@ void write_output_parallel(struct engine* e, const char* baseName,
hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
+ const size_t Nstars = e->s->nr_sparts;
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
int numFiles = 1;
struct part* parts = e->s->parts;
struct gpart* gparts = e->s->gparts;
struct gpart* dmparts = NULL;
+ struct spart* sparts = e->s->sparts;
static int outputCount = 0;
FILE* xmfFile = 0;
/* Number of unassociated gparts */
- const size_t Ndm = Ntot > 0 ? Ntot - Ngas : 0;
+ const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
/* File name */
char fileName[FILENAME_BUFFER_SIZE];
@@ -642,7 +674,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
/* Compute offset in the file and total number of
* particles */
- size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0};
+ size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0, 0, Nstars, 0};
long long N_total[NUM_PARTICLE_TYPES] = {0};
long long offset[NUM_PARTICLE_TYPES] = {0};
MPI_Exscan(&N, &offset, NUM_PARTICLE_TYPES, MPI_LONG_LONG_INT, MPI_SUM, comm);
@@ -816,11 +848,13 @@ void write_output_parallel(struct engine* e, const char* baseName,
/* Write DM particles */
Nparticles = Ndm;
darkmatter_write_particles(dmparts, list, &num_fields);
-
- /* Free temporary array */
- free(dmparts);
break;
+ case STAR:
+ Nparticles = Nstars;
+ star_write_particles(sparts, list, &num_fields);
+ break;
+
default:
error("Particle Type %d not yet supported. Aborting", ptype);
}
@@ -830,9 +864,12 @@ void write_output_parallel(struct engine* e, const char* baseName,
writeArray(e, h_grp, fileName, xmfFile, partTypeGroupName, list[i],
Nparticles, N_total[ptype], mpi_rank, offset[ptype],
internal_units, snapshot_units);
-
+
/* Free temporary array */
- free(dmparts);
+ if (dmparts) {
+ free(dmparts);
+ dmparts = 0;
+ }
/* Close particle group */
H5Gclose(h_grp);
diff --git a/src/parallel_io.h b/src/parallel_io.h
index e5b12aa50c30b4d63ccc81835d2d8454e01b3889..3f7395b8c17aaee8654d3b3f61e78d35cd27cb61 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -36,9 +36,11 @@
void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
- size_t* Ngas, size_t* Ngparts, int* periodic,
- int* flag_entropy, int mpi_rank, int mpi_size,
- MPI_Comm comm, MPI_Info info, int dry_run);
+ struct spart** sparts, size_t* Ngas, size_t* Ngparts,
+ size_t* Nsparts, int* periodic, int* flag_entropy,
+ int with_hydro, int with_gravity, int with_stars,
+ int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
+ int dry_run);
void write_output_parallel(struct engine* e, const char* baseName,
const struct UnitSystem* internal_units,