From 66d43308d1200354dd332f31548bc961edfee157 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 15 Mar 2016 14:54:00 +0100 Subject: [PATCH 01/17] Allow for more than 2^32 particles in single i/o mode --- src/single_io.c | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/single_io.c b/src/single_io.c index 59686a68b..c1c23fb07 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -204,7 +204,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, /* Make sure the chunks are not larger than the dataset */ if (chunk_shape[0] > N) chunk_shape[0] = N; - + /* Change shape of data space */ h_err = H5Sset_extent_simple(h_space, rank, shape, NULL); if (h_err < 0) { @@ -337,6 +337,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, double boxSize[3] = {0.0, -1.0, -1.0}; /* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/ int numParticles[NUM_PARTICLE_TYPES] = {0}; + int numParticles_highWord[NUM_PARTICLE_TYPES] = {0}; size_t Ndm; /* Open file */ @@ -365,9 +366,12 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, /* Read the relevant information and print status */ readAttribute(h_grp, "BoxSize", DOUBLE, boxSize); readAttribute(h_grp, "NumPart_Total", UINT, numParticles); + readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord); - *Ngas = numParticles[0]; - Ndm = numParticles[1]; + *Ngas = ((long long)numParticles[0]) + + ((long long)numParticles_highWord[0] << 32); + Ndm = ((long long)numParticles[1]) + + ((long long)numParticles_highWord[1] << 32); dim[0] = boxSize[0]; dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1]; dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2]; -- GitLab From 8423464356a07f2c41e721f0bf77a66fc6ccafbc Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 15 Mar 2016 14:58:01 +0100 Subject: [PATCH 02/17] Documentation fix --- src/single_io.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/single_io.c b/src/single_io.c index c1c23fb07..c5cbb5ac5 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -314,9 +314,9 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, * @param fileName The file to read. * @param dim (output) The dimension of the volume. * @param parts (output) Array of Gas particles. - * @param gparts (output) Array of DM particles. + * @param gparts (output) Array of #gpart particles. * @param Ngas (output) number of Gas particles read. - * @param Ngparts (output) The number of DM particles read. + * @param Ngparts (output) The number of #gpart read. * @param periodic (output) 1 if the volume is periodic, 0 if not. * * Opens the HDF5 file fileName and reads the particles contained -- GitLab From 1330a303aacf016d1e9d99f6bdafd565369d6af1 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Tue, 15 Mar 2016 15:18:45 +0100 Subject: [PATCH 03/17] First attempt at serial_io reads with multiple particle types --- src/serial_io.c | 117 +++++++++++++++++++++++++++++++++++------------- src/single_io.c | 19 ++++---- 2 files changed, 94 insertions(+), 42 deletions(-) diff --git a/src/serial_io.c b/src/serial_io.c index 8e63db5cf..8068d02db 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -396,9 +396,15 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, * * @param fileName The file to read. * @param dim (output) The dimension of the volume read from the file. - * @param parts (output) The array of #part read from the file. - * @param N (output) The number of particles read from the file. + * @param parts (output) The array of #part (gas particles) read from the file. + * @param gparts (output) The array of #gpart read from the file. + * @param Ngas (output) The number of #part read from the file on that node. + * @param Ngparts (output) The number of #gpart read from the file on that node. * @param periodic (output) 1 if the volume is periodic, 0 if not. + * @param mpi_rank The MPI rank of this node + * @param mpi_size The number of MPI ranks + * @param comm The MPI communicator + * @param info The MPI information object * * Opens the HDF5 file fileName and reads the particles contained * in the parts array. N is the returned number of particles found @@ -411,17 +417,18 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, * */ void read_ic_serial(char* fileName, double dim[3], struct part** parts, - size_t* N, int* periodic, int mpi_rank, int mpi_size, + struct gpart** gparts, size_t* Ngas, size_t* Ngparts, + int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info) { hid_t h_file = 0, h_grp = 0; - double boxSize[3] = {0.0, -1.0, -1.0}; /* GADGET has only cubic boxes (in cosmological mode) */ - int numParticles[6] = {0}; - /* GADGET has 6 particle types. We only keep the type 0*/ - int numParticles_highWord[6] = {0}; - long long offset = 0; - long long N_total = 0; - int rank; + double boxSize[3] = {0.0, -1.0, -1.0}; + /* GADGET has 6 particle types. We only keep the type 0 & 1 for now*/ + int numParticles[NUM_PARTICLE_TYPES] = {0}; + int numParticles_highWord[NUM_PARTICLE_TYPES] = {0}; + size_t N[NUM_PARTICLE_TYPES] = {0}; + long long N_total[NUM_PARTICLE_TYPES] = {0}; + long long offset[NUM_PARTICLE_TYPES] = {0}; /* First read some information about the content */ if (mpi_rank == 0) { @@ -453,8 +460,10 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, readAttribute(h_grp, "NumPart_Total", UINT, numParticles); readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord); - N_total = ((long long)numParticles[0]) + - ((long long)numParticles_highWord[0] << 32); + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + N_total[ptype] = ((long long)numParticles[ptype]) + + ((long long)numParticles_highWord[ptype] << 32); + dim[0] = boxSize[0]; dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1]; dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2]; @@ -474,42 +483,80 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, /* Now need to broadcast that information to all ranks. */ MPI_Bcast(periodic, 1, MPI_INT, 0, comm); - MPI_Bcast(&N_total, 1, MPI_LONG_LONG, 0, comm); + MPI_Bcast(&N_total, NUM_PARTICLE_TYPES, MPI_LONG_LONG, 0, comm); MPI_Bcast(dim, 3, MPI_DOUBLE, 0, comm); /* Divide the particles among the tasks. */ - offset = mpi_rank * N_total / mpi_size; - *N = (mpi_rank + 1) * N_total / mpi_size - offset; + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + offset[ptype] = mpi_rank * N_total[ptype] / mpi_size; + N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype]; + } - /* Allocate memory to store particles */ - if (posix_memalign((void*)parts, part_align, (*N) * sizeof(struct part)) != 0) + /* 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, *N * sizeof(struct part)); + bzero(*parts, *Ngas * sizeof(struct part)); + + /* Allocate memory to store all particles */ + 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) / */ /* (1024.*1024.)); */ + message("BoxSize = %lf", dim[0]); + message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); + /* Now loop over ranks and read the data */ - for (rank = 0; rank < mpi_size; ++rank) { + for (int rank = 0; rank < mpi_size; ++rank) { /* Is it this rank's turn to read ? */ if (rank == mpi_rank) { - + h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT); if (h_file < 0) error("Error while opening file '%s' on rank %d.", fileName, mpi_rank); - /* Open SPH particles group */ - /* message("Reading particle arrays..."); */ - h_grp = H5Gopen(h_file, "/PartType0", H5P_DEFAULT); - if (h_grp < 0) - error("Error while opening particle group on rank %d.\n", mpi_rank); - - /* Read particle fields into the particle structure */ - hydro_read_particles(h_grp, *N, N_total, offset, *parts); - - /* Close particle group */ - H5Gclose(h_grp); - + /* Loop over all particle types */ + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { + + /* Don't do anything if no particle of this kind */ + if (N[ptype] == 0) continue; + + /* Open the particle group in the file */ + char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; + snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", + ptype); + h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); + if (h_grp < 0) { + error("Error while opening particle group %s.", partTypeGroupName); + } + + + /* Read particle fields into the particle structure */ + switch (ptype) { + + case GAS: + hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype], *parts); + break; + + case DM: + darkmatter_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype], *gparts); + break; + + default: + error("Particle Type %d not yet supported. Aborting", ptype); + } + + /* Close particle group */ + H5Gclose(h_grp); + + } + /* Close file */ H5Fclose(h_file); } @@ -518,6 +565,12 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, MPI_Barrier(comm); } + /* Prepare the DM particles */ + prepare_dm_gparts(*gparts, Ndm); + + /* Now duplicate the hydro particle into gparts */ + duplicate_hydro_gparts(*parts, *gparts, *Ngas, Ndm); + /* message("Done Reading particles..."); */ } diff --git a/src/single_io.c b/src/single_io.c index c5cbb5ac5..32a4b374a 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -338,6 +338,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, /* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/ int numParticles[NUM_PARTICLE_TYPES] = {0}; int numParticles_highWord[NUM_PARTICLE_TYPES] = {0}; + size_t N[NUM_PARTICLE_TYPES] = {0}; size_t Ndm; /* Open file */ @@ -368,10 +369,10 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, readAttribute(h_grp, "NumPart_Total", UINT, numParticles); readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord); - *Ngas = ((long long)numParticles[0]) + - ((long long)numParticles_highWord[0] << 32); - Ndm = ((long long)numParticles[1]) + - ((long long)numParticles_highWord[1] << 32); + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + N[ptype] = ((long long)numParticles[ptype]) + + ((long long)numParticles_highWord[ptype] << 32); + dim[0] = boxSize[0]; dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1]; dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2]; @@ -382,16 +383,16 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, /* Close header */ H5Gclose(h_grp); - /* Total number of particles */ - *Ngparts = *Ngas + Ndm; - /* 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 SPH particles"); bzero(*parts, *Ngas * sizeof(struct part)); /* Allocate memory to store all particles */ + 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"); @@ -400,8 +401,6 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) / * (1024.*1024.)); */ - /* Open SPH particles group */ - /* message("Reading particle arrays..."); */ message("BoxSize = %lf", dim[0]); message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); @@ -409,7 +408,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { /* Don't do anything if no particle of this kind */ - if (numParticles[ptype] == 0) continue; + if (N[ptype] == 0) continue; /* Open the particle group in the file */ char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; -- GitLab From e48aa70dc05d3e21c684cc8fea8031db75e7be79 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 12:08:40 +0000 Subject: [PATCH 04/17] Implemented serial i/o reads --- examples/main.c | 2 +- src/common_io.h | 3 +++ src/serial_io.c | 5 ++++- src/serial_io.h | 3 ++- src/single_io.c | 3 --- 5 files changed, 10 insertions(+), 6 deletions(-) diff --git a/examples/main.c b/examples/main.c index c88f92a07..f747b6cb7 100644 --- a/examples/main.c +++ b/examples/main.c @@ -359,7 +359,7 @@ int main(int argc, char *argv[]) { read_ic_parallel(ICfileName, dim, &parts, &Ngas, &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL); #else - read_ic_serial(ICfileName, dim, &parts, &Ngas, &periodic, myrank, nr_nodes, + read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL); #endif #else diff --git a/src/common_io.h b/src/common_io.h index 2623a03f9..426fa6a01 100644 --- a/src/common_io.h +++ b/src/common_io.h @@ -70,6 +70,9 @@ enum PARTICLE_TYPE { NUM_PARTICLE_TYPES }; +#define FILENAME_BUFFER_SIZE 150 +#define PARTICLE_GROUP_BUFFER_SIZE 20 + hid_t hdf5Type(enum DATA_TYPE type); size_t sizeOfType(enum DATA_TYPE type); diff --git a/src/serial_io.c b/src/serial_io.c index 8068d02db..0d64db482 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -390,6 +390,9 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, /* Import the right hydro definition */ #include "hydro_io.h" +/* Import the right gravity definition */ +#include "gravity_io.h" + /** * @brief Reads an HDF5 initial condition file (GADGET-3 type) @@ -499,7 +502,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, bzero(*parts, *Ngas * sizeof(struct part)); /* Allocate memory to store all particles */ - Ndm = N[1]; + 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"); diff --git a/src/serial_io.h b/src/serial_io.h index 95f09f597..6cb48d26d 100644 --- a/src/serial_io.h +++ b/src/serial_io.h @@ -32,7 +32,8 @@ #if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5) void read_ic_serial(char* fileName, double dim[3], struct part** parts, - size_t* N, int* periodic, int mpi_rank, int mpi_size, + struct gpart** gparts, size_t* Ngas, size_t* Ngparts, + int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info); void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, diff --git a/src/single_io.c b/src/single_io.c index 32a4b374a..373d68517 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -39,9 +39,6 @@ #include "common_io.h" #include "error.h" -#define FILENAME_BUFFER_SIZE 150 -#define PARTICLE_GROUP_BUFFER_SIZE 20 - /*----------------------------------------------------------------------------- * Routines reading an IC file *-----------------------------------------------------------------------------*/ -- GitLab From eaa0b8f7438604c476d0063314008b212aefac63 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 12:47:07 +0000 Subject: [PATCH 05/17] Correct alignment for gparts --- src/gravity/Default/gravity_part.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h index 7ce7b8189..aaaab2055 100644 --- a/src/gravity/Default/gravity_part.h +++ b/src/gravity/Default/gravity_part.h @@ -50,4 +50,4 @@ struct gpart { struct part* part; }; -} __attribute__((aligned(part_align))); +} __attribute__((aligned(gpart_align))); -- GitLab From 0dab2629d6f2524064a049192c6bd7c3910d6fe2 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 13:59:26 +0000 Subject: [PATCH 06/17] Initialise engine_rank earlier --- examples/main.c | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/main.c b/examples/main.c index f747b6cb7..f8e918752 100644 --- a/examples/main.c +++ b/examples/main.c @@ -51,6 +51,8 @@ #define ENGINE_POLICY engine_policy_none #endif +extern int engine_rank; + /** * @brief Main routine that loads a few particles and generates some output. * @@ -120,6 +122,8 @@ int main(int argc, char *argv[]) { if (myrank == 0) message("MPI is up and running with %i node(s).", nr_nodes); fflush(stdout); + engine_rank = myrank; + /* Set a default grid so that grid[0]*grid[1]*grid[2] == nr_nodes. */ factor(nr_nodes, &initial_partition.grid[0], &initial_partition.grid[1]); factor(nr_nodes / initial_partition.grid[1], &initial_partition.grid[0], -- GitLab From 0d96b83352d1da0153f39d12921ff9c2dae5997c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 14:00:01 +0000 Subject: [PATCH 07/17] Start the clock before MPI --- examples/main.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/examples/main.c b/examples/main.c index f8e918752..dd112d28c 100644 --- a/examples/main.c +++ b/examples/main.c @@ -102,6 +102,9 @@ int main(int argc, char *argv[]) { /* Choke on FP-exceptions. */ // feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW ); + /* Initialize CPU frequency, this also starts time. */ + clocks_set_cpufreq(cpufreq); + #ifdef WITH_MPI /* Start by initializing MPI. */ int res, prov; @@ -132,9 +135,6 @@ int main(int argc, char *argv[]) { &initial_partition.grid[1], &initial_partition.grid[0]); #endif - /* Initialize CPU frequency, this also starts time. */ - clocks_set_cpufreq(cpufreq); - /* Greeting message */ if (myrank == 0) greetings(); -- GitLab From ed3f2c46a1fba9702fceaba54220cbcf5b8f06a3 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 16:07:38 +0000 Subject: [PATCH 08/17] Print correct number of particles --- examples/main.c | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/main.c b/examples/main.c index dd112d28c..b113621c7 100644 --- a/examples/main.c +++ b/examples/main.c @@ -380,6 +380,7 @@ int main(int argc, char *argv[]) { #if defined(WITH_MPI) long long N_long[2] = {Ngas, Ngpart}; MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + N_total[1] -= N_total[0]; if (myrank == 0) message("Read %lld gas particles and %lld DM particles from the ICs", N_total[0], N_total[1]); -- GitLab From 40df9f1d75db2b0b8bc793b85a9684e890c428fa Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 16:09:34 +0000 Subject: [PATCH 09/17] Make sure the right particle size is used --- src/common_io.c | 3 ++- src/single_io.c | 36 +++++++++++++++++++++--------------- 2 files changed, 23 insertions(+), 16 deletions(-) diff --git a/src/common_io.c b/src/common_io.c index f6a480333..4be334ccd 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -489,7 +489,8 @@ void prepare_dm_gparts(struct gpart* gparts, size_t Ndm) { for (size_t i = 0; i < Ndm; ++i) { /* 0 or negative ids are not allowed */ - if (gparts[i].id <= 0) error("0 or negative ID for DM particle"); + if (gparts[i].id <= 0) + error("0 or negative ID for DM particle %zd: ID=%lld", i, gparts[i].id); gparts[i].id = -gparts[i].id; } diff --git a/src/single_io.c b/src/single_io.c index 373d68517..f3082fd08 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -63,14 +63,14 @@ * Calls #error() if an error occurs. */ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, - int dim, char* part_c, enum DATA_IMPORTANCE importance) { + int dim, char* part_c, size_t partSize, + enum DATA_IMPORTANCE importance) { hid_t h_data = 0, h_err = 0, h_type = 0; htri_t exist = 0; void* temp; int i = 0; const size_t typeSize = sizeOfType(type); const size_t copySize = typeSize * dim; - const size_t partSize = sizeof(struct part); char* temp_c = 0; /* Check whether the dataspace exists or not */ @@ -155,14 +155,13 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, */ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, int N, int dim, char* part_c, - struct UnitSystem* us, + size_t partSize, struct UnitSystem* us, enum UnitConversionFactor convFactor) { hid_t h_data = 0, h_err = 0, h_space = 0, h_prop = 0; void* temp = 0; int i = 0, rank = 0; const size_t typeSize = sizeOfType(type); const size_t copySize = typeSize * dim; - const size_t partSize = sizeof(struct part); char* temp_c = 0; hsize_t shape[2]; hsize_t chunk_shape[2]; @@ -273,7 +272,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define readArray(grp, name, type, N, dim, part, N_total, offset, field, \ importance) \ readArrayBackEnd(grp, name, type, N, dim, (char*)(&(part[0]).field), \ - importance) + sizeof(part[0]),importance) /** * @brief A helper macro to call the readArrayBackEnd function more easily. @@ -298,7 +297,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, N_total, \ mpi_rank, offset, field, us, convFactor) \ writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, \ - (char*)(&(part[0]).field), us, convFactor) + (char*)(&(part[0]).field), sizeof(part[0]), us, convFactor) /* Import the right hydro definition */ #include "hydro_io.h" @@ -398,8 +397,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, /* 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("BoxSize = %lf", dim[0]); */ + /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */ /* Loop over all particle types */ for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { @@ -477,9 +476,8 @@ void write_output_single(struct engine* e, struct UnitSystem* us) { /* Number of particles of each type */ const size_t Ndm = Ntot - Ngas; - int numParticles[NUM_PARTICLE_TYPES] = /* Gadget-2 convention here */ - {Ngas, Ndm, 0}; /* Could use size_t instead */ - int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0}; + + long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0}; /* File name */ char fileName[FILENAME_BUFFER_SIZE]; @@ -521,19 +519,27 @@ void write_output_single(struct engine* e, struct UnitSystem* us) { /* Print the relevant information and print status */ writeAttribute(h_grp, "BoxSize", DOUBLE, e->s->dim, 3); - writeAttribute(h_grp, "NumPart_ThisFile", UINT, numParticles, - NUM_PARTICLE_TYPES); double dblTime = e->time; writeAttribute(h_grp, "Time", DOUBLE, &dblTime, 1); /* GADGET-2 legacy values */ + /* Number of particles of each type */ + unsigned int numParticles[NUM_PARTICLE_TYPES] = {0}; + unsigned int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0}; + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + numParticles[ptype] = (unsigned int) N_total[ptype]; + numParticlesHighWord[ptype] = (unsigned int) (N_total[ptype] >> 32); + } + writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total, + NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, NUM_PARTICLE_TYPES); - double MassTable[NUM_PARTICLE_TYPES] = {0., 0., 0., 0., 0., 0.}; + double MassTable[NUM_PARTICLE_TYPES] = {0}; writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES); - writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord, + unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0}; + writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1); -- GitLab From 85829ef097f78825ae2b8a9106a35caf914fac64 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 16:11:39 +0000 Subject: [PATCH 10/17] Write multi-types particles in serial i/o mode --- src/serial_io.c | 171 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 113 insertions(+), 58 deletions(-) diff --git a/src/serial_io.c b/src/serial_io.c index 0d64db482..05555f14a 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -68,7 +68,8 @@ */ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, long long N_total, long long offset, - char* part_c, enum DATA_IMPORTANCE importance) { + char* part_c, size_t partSize, + enum DATA_IMPORTANCE importance) { hid_t h_data = 0, h_err = 0, h_type = 0, h_memspace = 0, h_filespace = 0; hsize_t shape[2], offsets[2]; htri_t exist = 0; @@ -76,7 +77,6 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int i = 0, rank = 0; const size_t typeSize = sizeOfType(type); const size_t copySize = typeSize * dim; - const size_t partSize = sizeof(struct part); char* temp_c = 0; /* Check whether the dataspace exists or not */ @@ -269,7 +269,7 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name, */ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, int N, int dim, long long N_total, - int mpi_rank, long long offset, char* part_c, + int mpi_rank, long long offset, char* part_c, size_t partSize, struct UnitSystem* us, enum UnitConversionFactor convFactor) { @@ -279,7 +279,6 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, int i = 0, rank = 0; const size_t typeSize = sizeOfType(type); const size_t copySize = typeSize * dim; - const size_t partSize = sizeof(struct part); char* temp_c = 0; /* message("Writing '%s' array...", name); */ @@ -362,7 +361,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define readArray(grp, name, type, N, dim, part, N_total, offset, field, \ importance) \ readArrayBackEnd(grp, name, type, N, dim, N_total, offset, \ - (char*)(&(part[0]).field), importance) + (char*)(&(part[0]).field), sizeof(part[0]), importance) /** * @brief A helper macro to call the readArrayBackEnd function more easily. @@ -385,8 +384,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, N_total, \ mpi_rank, offset, field, us, convFactor) \ writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, N_total, \ - mpi_rank, offset, (char*)(&(part[0]).field), us, \ - convFactor) + mpi_rank, offset, (char*)(&(part[0]).field), \ + sizeof(part[0]), us, convFactor) /* Import the right hydro definition */ #include "hydro_io.h" @@ -511,8 +510,8 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, /* 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("BoxSize = %lf", dim[0]); */ + /* message("NumPart = [%zd, %zd] Total = %zd", *Ngas, Ndm, *Ngparts); */ /* Now loop over ranks and read the data */ for (int rank = 0; rank < mpi_size; ++rank) { @@ -581,7 +580,11 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, * @brief Writes an HDF5 output file (GADGET-3 type) with its XMF descriptor * * @param e The engine containing all the system. - * @param us The UnitSystem used for the conversion of units in the output + * @param us The UnitSystem used for the conversion of units in the output. + * @param mpi_rank The MPI rank of this node. + * @param mpi_size The number of MPI ranks. + * @param comm The MPI communicator. + * @param info The MPI information object * * Creates an HDF5 output file and writes the particles contained * in the engine. If such a file already exists, it is erased and replaced @@ -594,35 +597,33 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info) { hid_t h_file = 0, h_grp = 0, h_grpsph = 0; - int N = e->s->nr_parts; + const size_t Ngas = e->s->nr_parts; + const size_t Ntot = e->s->nr_gparts; + const size_t Ndm = Ntot - Ngas; int periodic = e->s->periodic; - int numParticles[6] = {N, 0}; - int numParticlesHighWord[6] = {0}; - unsigned int flagEntropy[6] = {0}; - long long N_total = 0, offset = 0; - double offset_d = 0., N_d = 0., N_total_d = 0.; int numFiles = 1; - int rank = 0; struct part* parts = e->s->parts; - FILE* xmfFile = 0; + struct gpart* gparts = e->s->gparts; + struct gpart* dmparts = NULL; static int outputCount = 0; + FILE* xmfFile = 0; /* File name */ - char fileName[200]; - sprintf(fileName, "output_%03i.hdf5", outputCount); + char fileName[FILENAME_BUFFER_SIZE]; + snprintf(fileName, FILENAME_BUFFER_SIZE, "output_%03i.hdf5", outputCount); /* Compute offset in the file and total number of particles */ - /* Done using double to allow for up to 2^50=10^15 particles */ - N_d = (double)N; - MPI_Exscan(&N_d, &offset_d, 1, MPI_DOUBLE, MPI_SUM, comm); - N_total_d = offset_d + N_d; - MPI_Bcast(&N_total_d, 1, MPI_DOUBLE, mpi_size - 1, comm); - if (N_total_d > 1.e15) - error( - "Error while computing the offset for parallel output: Simulation has " - "more than 10^15 particles.\n"); - N_total = (long long)N_total_d; - offset = (long long)offset_d; + size_t N[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 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, MPI_SUM, comm); + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + N_total[ptype] = offset[ptype] + N[ptype]; + + /* The last rank now has the correct N_total. Let's broadcast from there */ + MPI_Bcast(&N_total, 6, MPI_LONG_LONG, mpi_size - 1, comm); + + /* Now everybody konws its offset and the total number of particles of each type */ /* Do common stuff first */ if (mpi_rank == 0) { @@ -634,7 +635,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, xmfFile = prepareXMFfile(); /* Write the part corresponding to this specific output */ - writeXMFheader(xmfFile, N_total, fileName, e->time); + //writeXMFheader(xmfFile, N_total, fileName, e->time); /* Open file */ /* message("Opening file '%s'.", fileName); */ @@ -666,15 +667,20 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, writeAttribute(h_grp, "Time", DOUBLE, &dblTime, 1); /* GADGET-2 legacy values */ - numParticles[0] = (unsigned int)N_total; - writeAttribute(h_grp, "NumPart_ThisFile", UINT, numParticles, 6); - writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6); - numParticlesHighWord[0] = (unsigned int)(N_total >> 32); - writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, - 6); + /* Number of particles of each type */ + unsigned int numParticles[NUM_PARTICLE_TYPES] = {0}; + unsigned int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0}; + for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + numParticles[ptype] = (unsigned int) N_total[ptype]; + numParticlesHighWord[ptype] = (unsigned int) (N_total[ptype] >> 32); + } + writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total, NUM_PARTICLE_TYPES); + writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, NUM_PARTICLE_TYPES); + writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, NUM_PARTICLE_TYPES); double MassTable[6] = {0., 0., 0., 0., 0., 0.}; - writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6); - writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, 6); + writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES); + unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0}; + writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1); /* Close header */ @@ -692,21 +698,32 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, /* Print the system of Units */ writeUnitSystem(h_file, us); - /* Create SPH particles group */ - /* message("Writing particle arrays..."); */ - h_grp = - H5Gcreate(h_file, "/PartType0", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); - if (h_grp < 0) error("Error while creating particle group.\n"); + /* Loop over all particle types */ + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { - /* Close particle group */ - H5Gclose(h_grp); + /* Don't do anything if no particle of this kind */ + if (numParticles[ptype] == 0) continue; + /* Open the particle group in the file */ + char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; + snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", + ptype); + h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT, + H5P_DEFAULT); + if (h_grp < 0) { + error("Error while creating particle group.\n"); + } + + /* Close particle group */ + H5Gclose(h_grp); + } + /* Close file */ H5Fclose(h_file); } /* Now loop over ranks and write the data */ - for (rank = 0; rank < mpi_size; ++rank) { + for (int rank = 0; rank < mpi_size; ++rank) { /* Is it this rank's turn to write ? */ if (rank == mpi_rank) { @@ -715,19 +732,57 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, if (h_file < 0) error("Error while opening file '%s' on rank %d.", fileName, mpi_rank); - /* Open SPH particles group */ - /* message("Reading particle arrays..."); */ - h_grp = H5Gopen(h_file, "/PartType0", H5P_DEFAULT); - if (h_grp < 0) - error("Error while opening particle group on rank %d.\n", mpi_rank); + /* Loop over all particle types */ + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { - /* Write particle fields from the particle structure */ - hydro_write_particles(h_grp, fileName, xmfFile, N, N_total, mpi_rank, - offset, parts, us); + /* Don't do anything if no particle of this kind */ + if (N[ptype] == 0) continue; - /* Close particle group */ - H5Gclose(h_grp); + /* Open the particle group in the file */ + char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; + snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", + ptype); + h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); + if (h_grp < 0) { + error("Error while opening particle group %s.", partTypeGroupName); + } + + /* Read particle fields into the particle structure */ + switch (ptype) { + + + case GAS: + hydro_write_particles(h_grp, fileName, xmfFile, N[ptype], N_total[ptype], mpi_rank, + offset[ptype], parts, us); + + break; + + case DM: + /* Allocate temporary array */ + if (posix_memalign((void*)&dmparts, gpart_align, + Ndm * sizeof(struct gpart)) != 0) + error("Error while allocating temporart memory for DM particles"); + bzero(dmparts, Ndm * sizeof(struct gpart)); + + /* Collect the DM particles from gpart */ + collect_dm_gparts(gparts, Ntot, dmparts, Ndm); + + /* Write DM particles */ + darkmatter_write_particles(h_grp, fileName, xmfFile, N[ptype], N_total[ptype], mpi_rank, offset[ptype], dmparts, us); + + /* Free temporary array */ + free(dmparts); + break; + + default: + error("Particle Type %d not yet supported. Aborting", ptype); + } + + /* Close particle group */ + H5Gclose(h_grp); + } + /* Close file */ H5Fclose(h_file); } -- GitLab From 2e25a3f75b963dd33b05d3ec1ef6c012cd649ac1 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Wed, 16 Mar 2016 17:15:19 +0100 Subject: [PATCH 11/17] Code formatting --- examples/main.c | 14 ++- src/engine.c | 2 +- src/serial_io.c | 233 +++++++++++++++++++++++++----------------------- src/serial_io.h | 6 +- src/single_io.c | 17 ++-- 5 files changed, 138 insertions(+), 134 deletions(-) diff --git a/examples/main.c b/examples/main.c index b113621c7..898f093e9 100644 --- a/examples/main.c +++ b/examples/main.c @@ -51,8 +51,6 @@ #define ENGINE_POLICY engine_policy_none #endif -extern int engine_rank; - /** * @brief Main routine that loads a few particles and generates some output. * @@ -99,8 +97,8 @@ int main(int argc, char *argv[]) { #endif #endif -/* Choke on FP-exceptions. */ -// feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW ); + /* Choke on FP-exceptions. */ + // feenableexcept( FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW ); /* Initialize CPU frequency, this also starts time. */ clocks_set_cpufreq(cpufreq); @@ -125,8 +123,6 @@ int main(int argc, char *argv[]) { if (myrank == 0) message("MPI is up and running with %i node(s).", nr_nodes); fflush(stdout); - engine_rank = myrank; - /* Set a default grid so that grid[0]*grid[1]*grid[2] == nr_nodes. */ factor(nr_nodes, &initial_partition.grid[0], &initial_partition.grid[1]); factor(nr_nodes / initial_partition.grid[1], &initial_partition.grid[0], @@ -363,8 +359,8 @@ int main(int argc, char *argv[]) { read_ic_parallel(ICfileName, dim, &parts, &Ngas, &periodic, myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL); #else - read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic, myrank, nr_nodes, - MPI_COMM_WORLD, MPI_INFO_NULL); + read_ic_serial(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic, + myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL); #endif #else read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic); @@ -388,7 +384,7 @@ int main(int argc, char *argv[]) { N_total[0] = Ngas; N_total[1] = Ngpart - Ngas; message("Read %lld gas particles and %lld DM particles from the ICs", - N_total[0], N_total[1]); + N_total[0], N_total[1]); #endif /* Apply h scaling */ diff --git a/src/engine.c b/src/engine.c index b76585353..b7a15fd66 100644 --- a/src/engine.c +++ b/src/engine.c @@ -743,7 +743,7 @@ int engine_exchange_strays(struct engine *e, int offset, size_t *ind, *neighbours * * Here we construct all the tasks for all possible neighbouring non-empty - * local cells in the hierarchy. No dependencies are being added thus far. + * local cells in the hierarchy. No dependencies are being added thus far. * Additional loop over neighbours can later be added by simply duplicating * all the tasks created by this function. * diff --git a/src/serial_io.c b/src/serial_io.c index 05555f14a..99c927d83 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -69,7 +69,7 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, long long N_total, long long offset, char* part_c, size_t partSize, - enum DATA_IMPORTANCE importance) { + enum DATA_IMPORTANCE importance) { hid_t h_data = 0, h_err = 0, h_type = 0, h_memspace = 0, h_filespace = 0; hsize_t shape[2], offsets[2]; htri_t exist = 0; @@ -269,8 +269,8 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, char* name, */ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, enum DATA_TYPE type, int N, int dim, long long N_total, - int mpi_rank, long long offset, char* part_c, size_t partSize, - struct UnitSystem* us, + int mpi_rank, long long offset, char* part_c, + size_t partSize, struct UnitSystem* us, enum UnitConversionFactor convFactor) { hid_t h_data = 0, h_err = 0, h_memspace = 0, h_filespace = 0; @@ -361,7 +361,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define readArray(grp, name, type, N, dim, part, N_total, offset, field, \ importance) \ readArrayBackEnd(grp, name, type, N, dim, N_total, offset, \ - (char*)(&(part[0]).field), sizeof(part[0]), importance) + (char*)(&(part[0]).field), sizeof(part[0]), importance) /** * @brief A helper macro to call the readArrayBackEnd function more easily. @@ -384,15 +384,14 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, N_total, \ mpi_rank, offset, field, us, convFactor) \ writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, N_total, \ - mpi_rank, offset, (char*)(&(part[0]).field), \ - sizeof(part[0]), us, convFactor) + mpi_rank, offset, (char*)(&(part[0]).field), \ + sizeof(part[0]), us, convFactor) /* Import the right hydro definition */ #include "hydro_io.h" /* Import the right gravity definition */ #include "gravity_io.h" - /** * @brief Reads an HDF5 initial condition file (GADGET-3 type) * @@ -419,9 +418,9 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, * */ void read_ic_serial(char* fileName, double dim[3], struct part** parts, - struct gpart** gparts, size_t* Ngas, size_t* Ngparts, - int* periodic, int mpi_rank, int mpi_size, - MPI_Comm comm, MPI_Info info) { + struct gpart** gparts, size_t* Ngas, size_t* Ngparts, + int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm, + MPI_Info info) { 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}; @@ -462,9 +461,9 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, readAttribute(h_grp, "NumPart_Total", UINT, numParticles); readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord); - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) N_total[ptype] = ((long long)numParticles[ptype]) + - ((long long)numParticles_highWord[ptype] << 32); + ((long long)numParticles_highWord[ptype] << 32); dim[0] = boxSize[0]; dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1]; @@ -489,36 +488,38 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, MPI_Bcast(dim, 3, MPI_DOUBLE, 0, comm); /* Divide the particles among the tasks. */ - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { offset[ptype] = mpi_rank * N_total[ptype] / mpi_size; N[ptype] = (mpi_rank + 1) * N_total[ptype] / mpi_size - offset[ptype]; } /* Allocate memory to store SPH particles */ *Ngas = N[0]; - if (posix_memalign((void*)parts, part_align, (*Ngas) * sizeof(struct part)) != 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) + 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); */ - + /* Now loop over ranks and read the data */ for (int rank = 0; rank < mpi_size; ++rank) { /* Is it this rank's turn to read ? */ if (rank == mpi_rank) { - + h_file = H5Fopen(fileName, H5F_ACC_RDONLY, H5P_DEFAULT); if (h_file < 0) error("Error while opening file '%s' on rank %d.", fileName, mpi_rank); @@ -526,39 +527,39 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, /* Loop over all particle types */ for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { - /* Don't do anything if no particle of this kind */ - if (N[ptype] == 0) continue; - - /* Open the particle group in the file */ - char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; - snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", - ptype); - h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); - if (h_grp < 0) { - error("Error while opening particle group %s.", partTypeGroupName); - } - - - /* Read particle fields into the particle structure */ - switch (ptype) { - - case GAS: - hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype], *parts); - break; - - case DM: - darkmatter_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype], *gparts); - break; - - default: - error("Particle Type %d not yet supported. Aborting", ptype); - } - - /* Close particle group */ - H5Gclose(h_grp); - + /* Don't do anything if no particle of this kind */ + if (N[ptype] == 0) continue; + + /* Open the particle group in the file */ + char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; + snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", + ptype); + h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); + if (h_grp < 0) { + error("Error while opening particle group %s.", partTypeGroupName); + } + + /* Read particle fields into the particle structure */ + switch (ptype) { + + case GAS: + hydro_read_particles(h_grp, N[ptype], N_total[ptype], offset[ptype], + *parts); + break; + + case DM: + darkmatter_read_particles(h_grp, N[ptype], N_total[ptype], + offset[ptype], *gparts); + break; + + default: + error("Particle Type %d not yet supported. Aborting", ptype); + } + + /* Close particle group */ + H5Gclose(h_grp); } - + /* Close file */ H5Fclose(h_file); } @@ -617,13 +618,14 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, 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, MPI_SUM, comm); - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) N_total[ptype] = offset[ptype] + N[ptype]; /* The last rank now has the correct N_total. Let's broadcast from there */ MPI_Bcast(&N_total, 6, MPI_LONG_LONG, mpi_size - 1, comm); - /* Now everybody konws its offset and the total number of particles of each type */ + /* Now everybody konws its offset and the total number of particles of each + * type */ /* Do common stuff first */ if (mpi_rank == 0) { @@ -635,7 +637,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, xmfFile = prepareXMFfile(); /* Write the part corresponding to this specific output */ - //writeXMFheader(xmfFile, N_total, fileName, e->time); + // writeXMFheader(xmfFile, N_total, fileName, e->time); /* Open file */ /* message("Opening file '%s'.", fileName); */ @@ -670,17 +672,21 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, /* Number of particles of each type */ unsigned int numParticles[NUM_PARTICLE_TYPES] = {0}; unsigned int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0}; - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { - numParticles[ptype] = (unsigned int) N_total[ptype]; - numParticlesHighWord[ptype] = (unsigned int) (N_total[ptype] >> 32); - } - writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total, NUM_PARTICLE_TYPES); - writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, NUM_PARTICLE_TYPES); - writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, NUM_PARTICLE_TYPES); + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + numParticles[ptype] = (unsigned int)N_total[ptype]; + numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32); + } + writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total, + NUM_PARTICLE_TYPES); + writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, + NUM_PARTICLE_TYPES); + writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord, + NUM_PARTICLE_TYPES); double MassTable[6] = {0., 0., 0., 0., 0., 0.}; writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES); unsigned int flagEntropy[NUM_PARTICLE_TYPES] = {0}; - writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, NUM_PARTICLE_TYPES); + writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, + NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1); /* Close header */ @@ -707,17 +713,17 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, /* Open the particle group in the file */ char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", - ptype); + ptype); h_grp = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT, - H5P_DEFAULT); + H5P_DEFAULT); if (h_grp < 0) { - error("Error while creating particle group.\n"); + error("Error while creating particle group.\n"); } - + /* Close particle group */ H5Gclose(h_grp); } - + /* Close file */ H5Fclose(h_file); } @@ -735,54 +741,55 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, /* Loop over all particle types */ for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { - /* Don't do anything if no particle of this kind */ - if (N[ptype] == 0) continue; - - /* Open the particle group in the file */ - char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; - snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", - ptype); - h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); - if (h_grp < 0) { - error("Error while opening particle group %s.", partTypeGroupName); - } - - /* Read particle fields into the particle structure */ - switch (ptype) { - - - case GAS: - hydro_write_particles(h_grp, fileName, xmfFile, N[ptype], N_total[ptype], mpi_rank, - offset[ptype], parts, us); - - break; - - case DM: - /* Allocate temporary array */ - if (posix_memalign((void*)&dmparts, gpart_align, - Ndm * sizeof(struct gpart)) != 0) - error("Error while allocating temporart memory for DM particles"); - bzero(dmparts, Ndm * sizeof(struct gpart)); - - /* Collect the DM particles from gpart */ - collect_dm_gparts(gparts, Ntot, dmparts, Ndm); - - /* Write DM particles */ - darkmatter_write_particles(h_grp, fileName, xmfFile, N[ptype], N_total[ptype], mpi_rank, offset[ptype], dmparts, us); - - /* Free temporary array */ - free(dmparts); - break; - - default: - error("Particle Type %d not yet supported. Aborting", ptype); - - } - - /* Close particle group */ - H5Gclose(h_grp); + /* Don't do anything if no particle of this kind */ + if (N[ptype] == 0) continue; + + /* Open the particle group in the file */ + char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; + snprintf(partTypeGroupName, PARTICLE_GROUP_BUFFER_SIZE, "/PartType%d", + ptype); + h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT); + if (h_grp < 0) { + error("Error while opening particle group %s.", partTypeGroupName); + } + + /* Read particle fields into the particle structure */ + switch (ptype) { + + case GAS: + hydro_write_particles(h_grp, fileName, xmfFile, N[ptype], + N_total[ptype], mpi_rank, offset[ptype], + parts, us); + + break; + + case DM: + /* Allocate temporary array */ + if (posix_memalign((void*)&dmparts, gpart_align, + Ndm * sizeof(struct gpart)) != 0) + error("Error while allocating temporart memory for DM particles"); + bzero(dmparts, Ndm * sizeof(struct gpart)); + + /* Collect the DM particles from gpart */ + collect_dm_gparts(gparts, Ntot, dmparts, Ndm); + + /* Write DM particles */ + darkmatter_write_particles(h_grp, fileName, xmfFile, N[ptype], + N_total[ptype], mpi_rank, offset[ptype], + dmparts, us); + + /* Free temporary array */ + free(dmparts); + break; + + default: + error("Particle Type %d not yet supported. Aborting", ptype); + } + + /* Close particle group */ + H5Gclose(h_grp); } - + /* Close file */ H5Fclose(h_file); } diff --git a/src/serial_io.h b/src/serial_io.h index 6cb48d26d..5a34d420c 100644 --- a/src/serial_io.h +++ b/src/serial_io.h @@ -32,9 +32,9 @@ #if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5) void read_ic_serial(char* fileName, double dim[3], struct part** parts, - struct gpart** gparts, size_t* Ngas, size_t* Ngparts, - int* periodic, int mpi_rank, int mpi_size, - MPI_Comm comm, MPI_Info info); + struct gpart** gparts, size_t* Ngas, size_t* Ngparts, + int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm, + MPI_Info info); void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info); diff --git a/src/single_io.c b/src/single_io.c index f3082fd08..cb1ac496d 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -64,7 +64,7 @@ */ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N, int dim, char* part_c, size_t partSize, - enum DATA_IMPORTANCE importance) { + enum DATA_IMPORTANCE importance) { hid_t h_data = 0, h_err = 0, h_type = 0; htri_t exist = 0; void* temp; @@ -272,7 +272,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define readArray(grp, name, type, N, dim, part, N_total, offset, field, \ importance) \ readArrayBackEnd(grp, name, type, N, dim, (char*)(&(part[0]).field), \ - sizeof(part[0]),importance) + sizeof(part[0]), importance) /** * @brief A helper macro to call the readArrayBackEnd function more easily. @@ -297,7 +297,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name, #define writeArray(grp, fileName, xmfFile, name, type, N, dim, part, N_total, \ mpi_rank, offset, field, us, convFactor) \ writeArrayBackEnd(grp, fileName, xmfFile, name, type, N, dim, \ - (char*)(&(part[0]).field), sizeof(part[0]), us, convFactor) + (char*)(&(part[0]).field), sizeof(part[0]), us, \ + convFactor) /* Import the right hydro definition */ #include "hydro_io.h" @@ -365,7 +366,7 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, readAttribute(h_grp, "NumPart_Total", UINT, numParticles); readAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticles_highWord); - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) N[ptype] = ((long long)numParticles[ptype]) + ((long long)numParticles_highWord[ptype] << 32); @@ -526,10 +527,10 @@ void write_output_single(struct engine* e, struct UnitSystem* us) { /* Number of particles of each type */ unsigned int numParticles[NUM_PARTICLE_TYPES] = {0}; unsigned int numParticlesHighWord[NUM_PARTICLE_TYPES] = {0}; - for(int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { - numParticles[ptype] = (unsigned int) N_total[ptype]; - numParticlesHighWord[ptype] = (unsigned int) (N_total[ptype] >> 32); - } + for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ++ptype) { + numParticles[ptype] = (unsigned int)N_total[ptype]; + numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32); + } writeAttribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total, NUM_PARTICLE_TYPES); writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, -- GitLab From 93b372e5dab8f44ed7c0c61a1aeeaf3cd428457c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 17 Mar 2016 08:06:15 +0000 Subject: [PATCH 12/17] Uninitialised variables triggering warnings in DDT --- examples/main.c | 2 +- src/engine.c | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/main.c b/examples/main.c index b113621c7..1b52ecdda 100644 --- a/examples/main.c +++ b/examples/main.c @@ -81,7 +81,7 @@ int main(int argc, char *argv[]) { int nr_nodes = 1, myrank = 0; FILE *file_thread; int with_outputs = 1; - int verbose = 0, talking; + int verbose = 0, talking = 0; unsigned long long cpufreq = 0; #ifdef WITH_MPI diff --git a/src/engine.c b/src/engine.c index b76585353..63cde4239 100644 --- a/src/engine.c +++ b/src/engine.c @@ -1974,6 +1974,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads, e->dt_max = dt_max; e->file_stats = NULL; e->verbose = verbose; + e->count_step = 0; e->wallclock_time = 0.f; engine_rank = nodeID; -- GitLab From 3d5a42df0f4567fa738c1d2d3884b5c38f55721c Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 17 Mar 2016 08:14:37 +0000 Subject: [PATCH 13/17] Better check for whether or not to create the HDF5 groups for a given particle type --- src/serial_io.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/serial_io.c b/src/serial_io.c index 05555f14a..e4b92e5d0 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -702,7 +702,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { /* Don't do anything if no particle of this kind */ - if (numParticles[ptype] == 0) continue; + if (N_total[ptype] == 0) continue; /* Open the particle group in the file */ char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; @@ -736,7 +736,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { /* Don't do anything if no particle of this kind */ - if (N[ptype] == 0) continue; + if (N_total[ptype] == 0) continue; /* Open the particle group in the file */ char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; -- GitLab From 6bf681f3c8349e65fadee1e8ec4d65264e606908 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Thu, 17 Mar 2016 08:19:46 +0000 Subject: [PATCH 14/17] Better check for whether or not to create the HDF5 groups for a given particle type --- src/serial_io.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/serial_io.c b/src/serial_io.c index 9804a7012..c9b08c6e6 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -742,7 +742,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) { /* Don't do anything if no particle of this kind */ - if (N[ptype] == 0) continue; + if (N_total[ptype] == 0) continue; /* Open the particle group in the file */ char partTypeGroupName[PARTICLE_GROUP_BUFFER_SIZE]; -- GitLab From 669d099966dc5cc8a53c46124db33233ecfae3a3 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Fri, 18 Mar 2016 13:35:50 +0100 Subject: [PATCH 15/17] Temporary fix to preserve master --- examples/main.c | 20 ++++++++++++++++++++ src/serial_io.c | 8 +++++++- src/single_io.c | 6 +++++- 3 files changed, 32 insertions(+), 2 deletions(-) diff --git a/examples/main.c b/examples/main.c index 180505f60..9fb55e099 100644 --- a/examples/main.c +++ b/examples/main.c @@ -387,6 +387,26 @@ int main(int argc, char *argv[]) { N_total[0], N_total[1]); #endif + /* MATTHIEU: Temporary fix to preserve master */ + free(gparts); + Ngpart = 0; +#if defined(WITH_MPI) + N_long[0] = Ngas; + N_long[1] = Ngpart; + MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + if (myrank == 0) + message( + "AFTER FIX: Read %lld gas particles and %lld DM particles from the ICs", + N_total[0], N_total[1]); +#else + N_total[0] = Ngas; + N_total[1] = Ngpart; + message( + "AFTER FIX: Read %lld gas particles and %lld DM particles from the ICs", + N_total[0], N_total[1]); +#endif + /* MATTHIEU: End temporary fix */ + /* Apply h scaling */ if (scaling != 1.0) for (size_t k = 0; k < Ngas; k++) parts[k].h *= scaling; diff --git a/src/serial_io.c b/src/serial_io.c index c9b08c6e6..3fb1a1d1b 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -600,7 +600,6 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, hid_t h_file = 0, h_grp = 0, h_grpsph = 0; const size_t Ngas = e->s->nr_parts; const size_t Ntot = e->s->nr_gparts; - const size_t Ndm = Ntot - Ngas; int periodic = e->s->periodic; int numFiles = 1; struct part* parts = e->s->parts; @@ -609,6 +608,13 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank, static int outputCount = 0; FILE* xmfFile = 0; + /* Number of particles of each type */ + //const size_t Ndm = Ntot - Ngas; + + /* MATTHIEU: Temporary fix to preserve master */ + const size_t Ndm = Ntot > 0 ? Ntot - Ngas: 0; + /* MATTHIEU: End temporary fix */ + /* File name */ char fileName[FILENAME_BUFFER_SIZE]; snprintf(fileName, FILENAME_BUFFER_SIZE, "output_%03i.hdf5", outputCount); diff --git a/src/single_io.c b/src/single_io.c index cb1ac496d..54ace5c4e 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -476,7 +476,11 @@ void write_output_single(struct engine* e, struct UnitSystem* us) { static int outputCount = 0; /* Number of particles of each type */ - const size_t Ndm = Ntot - Ngas; + //const size_t Ndm = Ntot - Ngas; + + /* MATTHIEU: Temporary fix to preserve master */ + const size_t Ndm = Ntot > 0 ? Ntot - Ngas: 0; + /* MATTHIEU: End temporary fix */ long long N_total[NUM_PARTICLE_TYPES] = {Ngas, Ndm, 0}; -- GitLab From ae7871200981b8cb6f8c96fd457c0646affe3554 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Fri, 18 Mar 2016 13:45:51 +0100 Subject: [PATCH 16/17] Added one command-line option '-g' to switch on gravity (i.e. not de-allocate the gparts AND add the right policy to the engine() --- examples/main.c | 41 ++++++++++++++++++++++++++--------------- 1 file changed, 26 insertions(+), 15 deletions(-) diff --git a/examples/main.c b/examples/main.c index 9fb55e099..92c68e4ce 100644 --- a/examples/main.c +++ b/examples/main.c @@ -79,6 +79,8 @@ int main(int argc, char *argv[]) { int nr_nodes = 1, myrank = 0; FILE *file_thread; int with_outputs = 1; + int with_gravity = 0; + int engine_policies = 0; int verbose = 0, talking = 0; unsigned long long cpufreq = 0; @@ -156,7 +158,7 @@ int main(int argc, char *argv[]) { bzero(&s, sizeof(struct space)); /* Parse the options */ - while ((c = getopt(argc, argv, "a:c:d:e:f:h:m:oP:q:R:s:t:v:w:y:z:")) != -1) + while ((c = getopt(argc, argv, "a:c:d:e:f:gh:m:oP:q:R:s:t:v:w:y:z:")) != -1) switch (c) { case 'a': if (sscanf(optarg, "%lf", &scaling) != 1) @@ -185,6 +187,9 @@ int main(int argc, char *argv[]) { case 'f': if (!strcpy(ICfileName, optarg)) error("Error parsing IC file name."); break; + case 'g': + with_gravity = 1; + break; case 'h': if (sscanf(optarg, "%llu", &cpufreq) != 1) error("Error parsing CPU frequency."); @@ -388,23 +393,26 @@ int main(int argc, char *argv[]) { #endif /* MATTHIEU: Temporary fix to preserve master */ - free(gparts); - Ngpart = 0; + if (!with_gravity) { + free(gparts); + Ngpart = 0; #if defined(WITH_MPI) - N_long[0] = Ngas; - N_long[1] = Ngpart; - MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); - if (myrank == 0) + N_long[0] = Ngas; + N_long[1] = Ngpart; + MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); + if (myrank == 0) + message( + "AFTER FIX: Read %lld gas particles and %lld DM particles from the " + "ICs", + N_total[0], N_total[1]); +#else + N_total[0] = Ngas; + N_total[1] = Ngpart; message( "AFTER FIX: Read %lld gas particles and %lld DM particles from the ICs", N_total[0], N_total[1]); -#else - N_total[0] = Ngas; - N_total[1] = Ngpart; - message( - "AFTER FIX: Read %lld gas particles and %lld DM particles from the ICs", - N_total[0], N_total[1]); #endif + } /* MATTHIEU: End temporary fix */ /* Apply h scaling */ @@ -469,12 +477,15 @@ int main(int argc, char *argv[]) { message("nr of cells at depth %i is %i.", data[0], data[1]); } + /* Construct the engine policy */ + engine_policies = ENGINE_POLICY | engine_policy_steal | engine_policy_hydro; + if (with_gravity) engine_policies |= engine_policy_external_gravity; + /* Initialize the engine with this space. */ if (myrank == 0) clocks_gettime(&tic); if (myrank == 0) message("nr_nodes is %i.", nr_nodes); engine_init(&e, &s, dt_max, nr_threads, nr_queues, nr_nodes, myrank, - ENGINE_POLICY | engine_policy_steal | engine_policy_hydro, 0, - time_end, dt_min, dt_max, talking); + engine_policies, 0, time_end, dt_min, dt_max, talking); if (myrank == 0 && verbose) { clocks_gettime(&toc); message("engine_init took %.3f %s.", clocks_diff(&tic, &toc), -- GitLab From 1a819dee75434a2e6c7ee5a8939cb3a8a2a691cb Mon Sep 17 00:00:00 2001 From: Matthieu Schaller Date: Fri, 18 Mar 2016 14:00:13 +0100 Subject: [PATCH 17/17] Unlink the gparts --- examples/main.c | 3 +++ 1 file changed, 3 insertions(+) diff --git a/examples/main.c b/examples/main.c index 92c68e4ce..2b33a12eb 100644 --- a/examples/main.c +++ b/examples/main.c @@ -1,3 +1,4 @@ + /******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk), @@ -395,6 +396,8 @@ int main(int argc, char *argv[]) { /* MATTHIEU: Temporary fix to preserve master */ if (!with_gravity) { free(gparts); + for(size_t k = 0; k < Ngas; ++k) + parts[k].gpart = NULL; Ngpart = 0; #if defined(WITH_MPI) N_long[0] = Ngas; -- GitLab