diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index 0fef11db624924bf543e880ed7d856a750b98865..2600d79a4a70de3e69fb76b82e11f6a4b72a73d2 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -19,18 +19,19 @@
##############################################################################
import h5py
+import sys
from numpy import *
# Generates a swift IC file containing a cartesian distribution of particles
# at a constant density and pressure in a cubic box
# Parameters
-periodic= 1 # 1 For periodic box
+periodic= 1 # 1 For periodic box
boxSize = 1.
-L = 50 # Number of particles along one axis
-rho = 2. # Density
-P = 1. # Pressure
-gamma = 5./3. # Gas adiabatic index
+L = int(sys.argv[1]) # Number of particles along one axis
+rho = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
fileName = "uniformBox.hdf5"
@@ -39,34 +40,6 @@ numPart = L**3
mass = boxSize**3 * rho / numPart
internalEnergy = P / ((gamma - 1.)*rho)
-#Generate particles
-coords = zeros((numPart, 3))
-v = zeros((numPart, 3))
-m = zeros((numPart, 1))
-h = zeros((numPart, 1))
-u = zeros((numPart, 1))
-ids = zeros((numPart, 1), dtype='L')
-
-for i in range(L):
- for j in range(L):
- for k in range(L):
- index = i*L*L + j*L + k
- x = i * boxSize / L + boxSize / (2*L)
- y = j * boxSize / L + boxSize / (2*L)
- z = k * boxSize / L + boxSize / (2*L)
- coords[index,0] = x
- coords[index,1] = y
- coords[index,2] = z
- v[index,0] = 0.
- v[index,1] = 0.
- v[index,2] = 0.
- m[index] = mass
- h[index] = 2.251 * boxSize / L
- u[index] = internalEnergy
- ids[index] = index
-
-
-
#--------------------------------------------------
#File
@@ -90,17 +63,39 @@ grp.attrs["PeriodicBoundariesOn"] = periodic
#Particle group
grp = file.create_group("/PartType0")
-ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
-ds[()] = coords
+
+v = zeros((numPart, 3))
ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), 2.251 * boxSize / L)
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False, dtype='L').reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
ds[()] = ids
+x = ids % L;
+y = ((ids - x) / L) % L;
+z = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 3))
+coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
file.close()
diff --git a/src/parallel_io.c b/src/parallel_io.c
index b1386330aeb48d5005492813edbf7116f3cf99f3..8ffd59a591f4f035d912156be08c33ef90b5d8ee 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -110,9 +110,9 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
offsets[0] = offset;
offsets[1] = 0;
} else {
- rank = 1;
+ rank = 2;
shape[0] = N;
- shape[1] = 0;
+ shape[1] = 1;
offsets[0] = offset;
offsets[1] = 0;
}
@@ -238,13 +238,13 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];
+ /* message("Found %d particles in a %speriodic box of size [%f %f %f].", */
+ /* N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+
/* Divide the particles among the tasks. */
offset = mpi_rank * N_total / mpi_size;
*N = (mpi_rank + 1) * N_total / mpi_size - offset;
- /* message("Found %d particles in a %speriodic box of size [%f %f %f].", */
- /* *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
-
/* Close header */
H5Gclose(h_grp);
diff --git a/src/serial_io.c b/src/serial_io.c
index 93583f7f90b02bad81864c9e4130193306b99fde..19b895ca5be87aac8e39c1c1f81ad911e294d522 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -87,13 +87,15 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
if (importance == COMPULSORY) {
error("Compulsory data set '%s' not present in the file.", name);
} else {
- for (i = 0; i < N; ++i) memset(part_c + i * partSize, 0, copySize);
+ for (i = 0; i < N; ++i)
+ memset(part_c + i * partSize, 0, copySize);
return;
}
}
- /* message( "Reading %s '%s' array...", importance == COMPULSORY ?
- * "compulsory": "optional ", name); */
+ /* message( "Reading %s '%s' array...", importance == COMPULSORY ? */
+ /* "compulsory": "optional ", name); */
+ /* fflush(stdout); */
/* Open data space */
h_data = H5Dopen1(grp, name);
@@ -117,9 +119,9 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
offsets[0] = offset;
offsets[1] = 0;
} else {
- rank = 1;
+ rank = 2;
shape[0] = N;
- shape[1] = 0;
+ shape[1] = 1;
offsets[0] = offset;
offsets[1] = 0;
}
@@ -131,6 +133,20 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
h_filespace = H5Dget_space(h_data);
H5Sselect_hyperslab(h_filespace, H5S_SELECT_SET, offsets, NULL, shape, NULL);
+ /* int rank_memspace = H5Sget_simple_extent_ndims(h_memspace); */
+ /* int rank_filespace = H5Sget_simple_extent_ndims(h_filespace); */
+
+ /* message("Memspace rank: %d", rank_memspace); */
+ /* message("Filespace rank: %d", rank_filespace); */
+ /* fflush(stdout); */
+
+ /* hsize_t dims_memspace[2], max_dims_memspace[2]; */
+ /* hsize_t dims_filespace[2], max_dims_filespace[2]; */
+
+ /* H5Sget_simple_extent_dims(h_memspace, dims_memspace, max_dims_memspace); */
+ /* H5Sget_simple_extent_dims(h_filespace, dims_filespace, max_dims_filespace);
+ */
+
/* Read HDF5 dataspace in temporary buffer */
/* Dirty version that happens to work for vectors but should be improved */
/* Using HDF5 dataspaces would be better */
@@ -196,11 +212,11 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
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};
+ 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;
@@ -235,14 +251,17 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
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);
+ N_total = ((long long) numParticles[0]) +
+ ((long long) numParticles_highWord[0] << 32);
dim[0] = boxSize[0];
dim[1] = (boxSize[1] < 0) ? boxSize[0] : boxSize[1];
dim[2] = (boxSize[2] < 0) ? boxSize[0] : boxSize[2];
- /* message("Found %d particles in a %speriodic box of size [%f %f %f].", */
- /* *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+ /* message("Found %lld particles in a %speriodic box of size [%f %f %f].",
+ */
+ /* N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
+
+ fflush(stdout);
/* Close header */
H5Gclose(h_grp);
@@ -264,8 +283,8 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
if (posix_memalign((void*)parts, part_align, (*N) * sizeof(struct part)) != 0)
error("Error while allocating memory for particles");
bzero(*parts, *N * sizeof(struct part));
- /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) /
- * (1024.*1024.)); */
+ /* message("Allocated %8.2f MB for particles.", *N * sizeof(struct part) / */
+ /* (1024.*1024.)); */
/* Now loop over ranks and read the data */
for (rank = 0; rank < mpi_size; ++rank) {
@@ -499,9 +518,9 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
hid_t h_file = 0, h_grp = 0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
- int numParticles[6] = {N, 0};
- int numParticlesHighWord[6] = {0};
- unsigned int flagEntropy[6] = {0};
+ 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;
@@ -516,7 +535,7 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
/* 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;
+ 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);
@@ -524,8 +543,8 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
error(
"Error while computing the offest for parallel output: Simulation has "
"more than 10^15 particles.\n");
- N_total = (long long)N_total_d;
- offset = (long long)offset_d;
+ N_total = (long long) N_total_d;
+ offset = (long long) offset_d;
/* Do common stuff first */
if (mpi_rank == 0) {
@@ -568,13 +587,13 @@ 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;
+ 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);
- double MassTable[6] = {0., 0., 0., 0., 0., 0.};
+ 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, "NumFilesPerSnapshot", INT, &numFiles, 1);
diff --git a/src/single_io.c b/src/single_io.c
index ddfb330370a1278f11f6696e81fe0066bd7b2562..0bdd9049b2dbb1d05479914e3f5e6791f5b4583d 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -83,7 +83,8 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
/* message("Optional data set '%s' not present. Zeroing this particle
* field...", name); */
- for (i = 0; i < N; ++i) memset(part_c + i * partSize, 0, copySize);
+ for (i = 0; i < N; ++i)
+ memset(part_c + i * partSize, 0, copySize);
return;
}
@@ -166,10 +167,10 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
int* periodic) {
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*/
+ 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*/
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -384,8 +385,8 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
hid_t h_file = 0, h_grp = 0;
int N = e->s->nr_parts;
int periodic = e->s->periodic;
- int numParticles[6] = {N, 0};
- int numParticlesHighWord[6] = {0};
+ int numParticles[6] = { N, 0 };
+ int numParticlesHighWord[6] = { 0 };
int numFiles = 1;
struct part* parts = e->s->parts;
FILE* xmfFile = 0;
@@ -437,7 +438,7 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6);
writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord,
6);
- double MassTable[6] = {0., 0., 0., 0., 0., 0.};
+ double MassTable[6] = { 0., 0., 0., 0., 0., 0. };
writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, 6);
writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord, 6);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);