diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index 8b668289525725bae3ca613feb76b331e7f1909f..0994e1c95e053defe7766122c52bc405c7239776 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -64,7 +64,7 @@ v = inputFile["/PartType0/Velocities"][:,:]
m = inputFile["/PartType0/Masses"][:]
h = inputFile["/PartType0/SmoothingLength"][:]
u = inputFile["/PartType0/InternalEnergy"][:]
-ids = np.array(range(np.size(u)), dtype='L')
+ids = np.array(range(np.size(u)), dtype='L') + 1
# Downsample
print "Downsampling..."
diff --git a/examples/BigPerturbedBox/makeIC_fcc.py b/examples/BigPerturbedBox/makeIC_fcc.py
index dfe2a41b7d8a17f3f7424294e11d944ae5ff2b94..13809b41b65890c51abfbd6db48640fb07b21623 100644
--- a/examples/BigPerturbedBox/makeIC_fcc.py
+++ b/examples/BigPerturbedBox/makeIC_fcc.py
@@ -111,6 +111,6 @@ ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids
+ds[()] = ids + 1
file.close()
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 533199869ee325dd75831410f3e9f1181ae30ebd..6aceeed559324f97b0b1e388ff0c3524498b52e4 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -82,7 +82,7 @@ for i in range(L):
else:
P = P + 3. + 4.*log(2.)
u[index] = P / ((gamma - 1.)*rho)
- ids[index] = partId
+ ids[index] = partId + 1
partId = partId + 1
diff --git a/examples/Makefile.am b/examples/Makefile.am
index b2c09e54800c556ba03bdb1e2436af3144263dcb..200d340be168651e6dfc91f4b89a0a343d34d9da 100644
--- a/examples/Makefile.am
+++ b/examples/Makefile.am
@@ -63,13 +63,15 @@ swift_fixdt_mpi_LDADD = ../src/.libs/libswiftsim_mpi.a $(HDF5_LDFLAGS) $(HDF5_L
# Scripts to generate ICs
EXTRA_DIST = UniformBox/makeIC.py \
- PerturbedBox/makeIC.py \
+ UniformDMBox/makeIC.py \
+ PerturbedBox/makeIC.py \
SedovBlast/makeIC.py SedovBlast/makeIC_fcc.py SedovBlast/solution.py \
SodShock/makeIC.py SodShock/solution.py SodShock/glass_001.hdf5 SodShock/glass_002.hdf5 SodShock/rhox.py \
CosmoVolume/getIC.sh \
BigCosmoVolume/makeIC.py \
BigPerturbedBox/makeIC_fcc.py \
- GreshoVortex/makeIC.py GreshoVortex/solution.py
+ GreshoVortex/makeIC.py GreshoVortex/solution.py \
+ MultiTypes/makeIC.py
# Scripts to plot task graphs
diff --git a/examples/MultiTypes/makeIC.py b/examples/MultiTypes/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..3a41910c22c260086b5384b248a5c86ab6340a5e
--- /dev/null
+++ b/examples/MultiTypes/makeIC.py
@@ -0,0 +1,136 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+import sys
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of DM particles
+# with a density of 1
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1.
+Lgas = int(sys.argv[1]) # Number of particles along one axis
+rhoGas = 2. # Density
+P = 1. # Pressure
+gamma = 5./3. # Gas adiabatic index
+rhoDM = 1.
+Ldm = int(sys.argv[2]) # Number of particles along one axis
+
+fileName = "multiTypes.hdf5"
+
+#---------------------------------------------------
+numGas = Lgas**3
+massGas = boxSize**3 * rhoGas / numGas
+internalEnergy = P / ((gamma - 1.)*rhoGas)
+
+numDM = Ldm**3
+massDM = boxSize**3 * rhoDM / numDM
+
+#--------------------------------------------------
+
+#File
+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_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numGas, numDM, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+
+# Gas Particle group
+grp = file.create_group("/PartType0")
+
+v = zeros((numGas, 3))
+ds = grp.create_dataset('Velocities', (numGas, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numGas, 1), massGas)
+ds = grp.create_dataset('Masses', (numGas,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+h = full((numGas, 1), 1.1255 * boxSize / Lgas)
+ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
+ds[()] = h
+h = zeros(1)
+
+u = full((numGas, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numGas,1), 'f')
+ds[()] = u
+u = zeros(1)
+
+ids = linspace(0, numGas, numGas, endpoint=False).reshape((numGas,1))
+ds = grp.create_dataset('ParticleIDs', (numGas, 1), 'L')
+ds[()] = ids + 1
+x = ids % Lgas;
+y = ((ids - x) / Lgas) % Lgas;
+z = (ids - x - Lgas * y) / Lgas**2;
+coords = zeros((numGas, 3))
+coords[:,0] = z[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+coords[:,1] = y[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+coords[:,2] = x[:,0] * boxSize / Lgas + boxSize / (2*Lgas)
+ds = grp.create_dataset('Coordinates', (numGas, 3), 'd')
+ds[()] = coords
+
+
+
+
+
+# DM Particle group
+grp = file.create_group("/PartType1")
+
+v = zeros((numDM, 3))
+ds = grp.create_dataset('Velocities', (numDM, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numDM, 1), massDM)
+ds = grp.create_dataset('Masses', (numDM,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+ids = linspace(0, numDM, numDM, endpoint=False).reshape((numDM,1))
+ds = grp.create_dataset('ParticleIDs', (numDM, 1), 'L')
+ds[()] = ids + Lgas**3 + 1
+x = ids % Ldm;
+y = ((ids - x) / Ldm) % Ldm;
+z = (ids - x - Ldm * y) / Ldm**2;
+coords = zeros((numDM, 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', (numDM, 3), 'd')
+ds[()] = coords
+
+file.close()
diff --git a/examples/PerturbedBox/makeIC.py b/examples/PerturbedBox/makeIC.py
index a5e831eca02463d287ce2c7748eb780ef66aeb33..69c1a69199c9a5262f5ae6c4e95ca14699300fd4 100644
--- a/examples/PerturbedBox/makeIC.py
+++ b/examples/PerturbedBox/makeIC.py
@@ -103,6 +103,6 @@ ds[()] = h
ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
ds[()] = u
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids
+ds[()] = ids + 1
file.close()
diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast/makeIC.py
index 04ff73a3cdc4d81cd6edb45ef66e4eb8b0df46a2..75ff81165df51780848e3d8ac679a6dbeb17a039 100644
--- a/examples/SedovBlast/makeIC.py
+++ b/examples/SedovBlast/makeIC.py
@@ -69,7 +69,7 @@ for i in range(L):
m[index] = mass
h[index] = 1.1255 * boxSize / L
u[index] = internalEnergy
- ids[index] = index
+ ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
u[index] = u[index] + E0 / (33. * mass)
print "Particle " , index , " set to detonate."
diff --git a/examples/SedovBlast/makeIC_fcc.py b/examples/SedovBlast/makeIC_fcc.py
index edff869b38d8077944574caf8273bca2b0ab6a33..17f07440909cb5478d09a5b7a1444c72af2f3a47 100644
--- a/examples/SedovBlast/makeIC_fcc.py
+++ b/examples/SedovBlast/makeIC_fcc.py
@@ -72,7 +72,7 @@ for i in range(L):
m[index] = mass
h[index] = 1.1255 * hbox
u[index] = internalEnergy
- ids[index] = index
+ ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
u[index] = u[index] + E0 / (28. * mass)
print "Particle " , index , " set to detonate."
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index 5bee473f06543d50412ccc476960504c10d56455..0ac0564116f8a6ceb57b4f41d23eb9907df0440d 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -86,8 +86,8 @@ h = append(h1, h2,0)
u = append(u1, u2,0)
m = append(m1, m2,0)
ids = zeros(numPart, dtype='L')
-for i in range(numPart):
- ids[i] = i
+for i in range(1, numPart+1):
+ ids[i-1] = i
#Final operations
h /= 2
diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index ea094e5a531456788dcbc35c098dc04e4d858da3..c175349e658799cbcb30dfe2619a1594bafc18b9 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -86,7 +86,7 @@ u = zeros(1)
ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
-ds[()] = ids
+ds[()] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
diff --git a/examples/UniformBox/makeICbig.py b/examples/UniformBox/makeICbig.py
index 28b4667cccf381532900c61228c24b3b6b1d62c1..e475fdcbd9f3c4811e3dcfdf20bbd321be3d8b29 100644
--- a/examples/UniformBox/makeICbig.py
+++ b/examples/UniformBox/makeICbig.py
@@ -98,7 +98,7 @@ for n in range(n_iterations):
u = zeros(1)
ids = linspace(offset, offset+N, N, endpoint=False).reshape((N,1))
- ds_id[offset:offset+N] = ids
+ ds_id[offset:offset+N] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
@@ -130,10 +130,8 @@ u = full((remainder, 1), internalEnergy)
ds_u[offset:offset+remainder] = u
u = zeros(1)
-print "Done", offset+remainder,"/", numPart
-
ids = linspace(offset, offset+remainder, remainder, endpoint=False).reshape((remainder,1))
-ds_id[offset:offset+remainder] = ids
+ds_id[offset:offset+remainder] = ids + 1
x = ids % L;
y = ((ids - x) / L) % L;
z = (ids - x - L * y) / L**2;
@@ -143,6 +141,7 @@ coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
ds_x[offset:offset+remainder,:] = coords
+print "Done", offset+remainder,"/", numPart
diff --git a/examples/UniformDMBox/makeIC.py b/examples/UniformDMBox/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..061b4d0ad1959d9e25356aff80e78adb9c1c4faa
--- /dev/null
+++ b/examples/UniformDMBox/makeIC.py
@@ -0,0 +1,86 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ #
+ # This program is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU Lesser General Public License as published
+ # by the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # This program is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU Lesser General Public License
+ # along with this program. If not, see <http://www.gnu.org/licenses/>.
+ #
+ ##############################################################################
+
+import h5py
+import sys
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of DM particles
+# with a density of 1
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1.
+rho = 1.
+L = int(sys.argv[1]) # Number of particles along one axis
+fileName = "uniformDMBox_%d.hdf5"%L
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, mass, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Particle group
+grp = file.create_group("/PartType1")
+
+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)
+
+ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+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/examples/main.c b/examples/main.c
index 0fce7500af2e3d4fe98feb7e4d3f484658bd2ded..c88f92a07a747c327692b5e0fbbc7dc07b93ac0c 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -58,8 +58,9 @@
int main(int argc, char *argv[]) {
- int c, icount, j, k, N, periodic = 1;
- long long N_total = -1;
+ int c, icount, periodic = 1;
+ size_t Ngas = 0, Ngpart = 0;
+ long long N_total[2] = {0, 0};
int nr_threads = 1, nr_queues = -1;
int dump_tasks = 0;
int data[2];
@@ -67,6 +68,7 @@ int main(int argc, char *argv[]) {
double h_max = -1.0, scaling = 1.0;
double time_end = DBL_MAX;
struct part *parts = NULL;
+ struct gpart *gparts = NULL;
struct space s;
struct engine e;
struct UnitSystem us;
@@ -354,14 +356,14 @@ 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, dim, &parts, &N, &periodic, myrank, nr_nodes,
+ read_ic_parallel(ICfileName, dim, &parts, &Ngas, &periodic, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL);
#else
- read_ic_serial(ICfileName, dim, &parts, &N, &periodic, myrank, nr_nodes,
+ read_ic_serial(ICfileName, dim, &parts, &Ngas, &periodic, myrank, nr_nodes,
MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
#else
- read_ic_single(ICfileName, dim, &parts, &N, &periodic);
+ read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic);
#endif
if (myrank == 0) {
@@ -372,24 +374,35 @@ int main(int argc, char *argv[]) {
}
#if defined(WITH_MPI)
- long long N_long = N;
- MPI_Reduce(&N_long, &N_total, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+ long long N_long[2] = {Ngas, Ngpart};
+ MPI_Reduce(&N_long, &N_total, 2, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+ if (myrank == 0)
+ message("Read %lld gas particles and %lld DM particles from the ICs",
+ N_total[0], N_total[1]);
#else
- N_total = N;
+ 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]);
#endif
- if (myrank == 0) message("Read %lld particles from the ICs", N_total);
/* Apply h scaling */
if (scaling != 1.0)
- for (k = 0; k < N; k++) parts[k].h *= scaling;
+ for (size_t k = 0; k < Ngas; k++) parts[k].h *= scaling;
/* Apply shift */
- if (shift[0] != 0 || shift[1] != 0 || shift[2] != 0)
- for (k = 0; k < N; k++) {
+ if (shift[0] != 0 || shift[1] != 0 || shift[2] != 0) {
+ for (size_t k = 0; k < Ngas; k++) {
parts[k].x[0] += shift[0];
parts[k].x[1] += shift[1];
parts[k].x[2] += shift[2];
}
+ for (size_t k = 0; k < Ngpart; k++) {
+ gparts[k].x[0] += shift[0];
+ gparts[k].x[1] += shift[1];
+ gparts[k].x[2] += shift[2];
+ }
+ }
/* Set default number of queues. */
if (nr_queues < 0) nr_queues = nr_threads;
@@ -399,7 +412,8 @@ int main(int argc, char *argv[]) {
/* Initialize the space with this data. */
if (myrank == 0) clocks_gettime(&tic);
- space_init(&s, dim, parts, N, periodic, h_max, myrank == 0);
+ space_init(&s, dim, parts, gparts, Ngas, Ngpart, periodic, h_max,
+ myrank == 0);
if (myrank == 0 && verbose) {
clocks_gettime(&toc);
message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -415,6 +429,7 @@ int main(int argc, char *argv[]) {
message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
s.cdim[1], s.cdim[2]);
message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
+ message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
message("maximum depth is %d.", s.maxdepth);
// message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
}
@@ -482,9 +497,10 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
message(
- "Running on %lld particles until t=%.3e with %i threads and %i "
- "queues (dt_min=%.3e, dt_max=%.3e)...",
- N_total, time_end, e.nr_threads, e.sched.nr_queues, e.dt_min, e.dt_max);
+ "Running on %lld gas particles and %lld DM particles until t=%.3e with "
+ "%i threads and %i queues (dt_min=%.3e, dt_max=%.3e)...",
+ N_total[0], N_total[1], time_end, e.nr_threads, e.sched.nr_queues,
+ e.dt_min, e.dt_max);
fflush(stdout);
}
@@ -499,7 +515,7 @@ int main(int argc, char *argv[]) {
clocks_getunit());
/* Let loose a runner on the space. */
- for (j = 0; !engine_is_done(&e); j++) {
+ for (int j = 0; !engine_is_done(&e); j++) {
/* Repartition the space amongst the nodes? */
#ifdef WITH_MPI
diff --git a/src/Makefile.am b/src/Makefile.am
index 8a4e7cdc4b849e8f99bdf9a3dbda9b30fa986332..f44d47819672d10445fd969fe2ff20dbcb49463b 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,9 +46,9 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel.h vector.h \
runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
- gravity.h \
- gravity/Default/gravity.h gravity/Default/runner_iact_grav.h \
- gravity/Default/gravity_part.h \
+ gravity.h gravity_io.h \
+ gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
+ gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
diff --git a/src/cell.c b/src/cell.c
index 3242c6b0b8226719d3f5311534f7945bb6c20a70..df11782048dfa80c697f53feefe8fabc104eb23b 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -345,9 +345,9 @@ void cell_split(struct cell *c) {
int i, j;
const int count = c->count, gcount = c->gcount;
- struct part* parts = c->parts;
- struct xpart* xparts = c->xparts;
- struct gpart* gparts = c->gparts;
+ struct part *parts = c->parts;
+ struct xpart *xparts = c->xparts;
+ struct gpart *gparts = c->gparts;
int left[8], right[8];
double pivot[3];
diff --git a/src/common_io.c b/src/common_io.c
index 0e9c5b6cd51ec64e04995cbc47be64b1f36bc131..f6a4803333581b69671e3adc223b46122ec5364c 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -359,7 +359,9 @@ FILE* prepareXMFfile() {
if (tempFile == NULL) error("Unable to open temporary file.");
- for (int i = 0; fgets(buffer, 1024, tempFile) != NULL && i < counter - 3; i++) {
+ int i = 0;
+ while (fgets(buffer, 1024, tempFile) != NULL && i < counter - 3) {
+ i++;
fprintf(xmfFile, "%s", buffer);
}
fprintf(xmfFile, "\n");
@@ -471,4 +473,90 @@ void writeXMFline(FILE* xmfFile, char* fileName, char* name, long long N,
fprintf(xmfFile, "</Attribute>\n");
}
+/**
+ * @brief Prepare the DM particles (in gparts) read in for the addition of the
+ * other particle types
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array
+ *
+ * @param gparts The array of #gpart freshly read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void prepare_dm_gparts(struct gpart* gparts, size_t Ndm) {
+
+ /* Let's give all these gparts a negative id */
+ 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");
+
+ gparts[i].id = -gparts[i].id;
+ }
+}
+
+/**
+ * @brief Copy every #part into the corresponding #gpart and link them.
+ *
+ * This function assumes that the DM particles are all at the start of the
+ * gparts array and adds the hydro particles afterwards
+ *
+ * @param parts The array of #part freshly read in.
+ * @param gparts The array of #gpart freshly read in with all the DM particles
+ *at the start
+ * @param Ngas The number of gas particles read in.
+ * @param Ndm The number of DM particles read in.
+ */
+void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
+ size_t Ngas, size_t Ndm) {
+
+ for (size_t i = 0; i < Ngas; ++i) {
+
+ /* Duplicate the crucial information */
+ gparts[i + Ndm].x[0] = parts[i].x[0];
+ gparts[i + Ndm].x[1] = parts[i].x[1];
+ gparts[i + Ndm].x[2] = parts[i].x[2];
+
+ gparts[i + Ndm].v_full[0] = parts[i].v[0];
+ gparts[i + Ndm].v_full[1] = parts[i].v[1];
+ gparts[i + Ndm].v_full[2] = parts[i].v[2];
+
+ gparts[i + Ndm].mass = parts[i].mass;
+
+ /* Link the particles */
+ gparts[i + Ndm].part = &parts[i];
+ parts[i].gpart = &gparts[i + Ndm];
+ }
+}
+
+/**
+ * @brief Copy every DM #gpart into the dmparts array.
+ *
+ * @param gparts The array of #gpart containing all particles.
+ * @param Ntot The number of #gpart.
+ * @param dmparts The array of #gpart containg DM particles to be filled.
+ * @param Ndm The number of DM particles.
+ */
+void collect_dm_gparts(struct gpart* gparts, size_t Ntot, struct gpart* dmparts,
+ size_t Ndm) {
+
+ size_t count = 0;
+
+ /* Loop over all gparts */
+ for (size_t i = 0; i < Ntot; ++i) {
+
+ /* And collect the DM ones */
+ if (gparts[i].id < 0) {
+ memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
+ dmparts[count].id = -dmparts[count].id;
+ count++;
+ }
+ }
+
+ /* Check that everything is fine */
+ if (count != Ndm)
+ error("Collected the wrong number of dm particles (%zd vs. %zd expected)",
+ count, Ndm);
+}
+
#endif
diff --git a/src/common_io.h b/src/common_io.h
index c221ad3aaf84deb83c0067ee4ece729b98003430..2623a03f9a25ce0e650dde4f698da6eb49177e26 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -24,6 +24,7 @@
#include "../config.h"
/* Includes. */
+#include "part.h"
#include "units.h"
#if defined(HAVE_HDF5)
@@ -55,9 +56,29 @@ enum DATA_IMPORTANCE {
OPTIONAL = 0
};
+/**
+ * @brief The different particle types present in a GADGET IC file
+ *
+ */
+enum PARTICLE_TYPE {
+ GAS = 0,
+ DM = 1,
+ BOUNDARY = 2,
+ DUMMY = 3,
+ STAR = 4,
+ BH = 5,
+ NUM_PARTICLE_TYPES
+};
+
hid_t hdf5Type(enum DATA_TYPE type);
size_t sizeOfType(enum DATA_TYPE type);
+void collect_dm_gparts(struct gpart* gparts, size_t Ntot, struct gpart* dmparts,
+ size_t Ndm);
+void prepare_dm_gparts(struct gpart* gparts, size_t Ndm);
+void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts,
+ size_t Ngas, size_t Ndm);
+
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
diff --git a/src/debug.c b/src/debug.c
index 48a0995a5d2992c2b23c7294833309508ee5b59d..4c1434118c98aab7def28d3a53493767d249d774 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -38,6 +38,8 @@
#error "Invalid choice of SPH variant"
#endif
+#include "./gravity/Default/gravity_debug.h"
+
/**
* @brief Looks for the particle with the given id and prints its information to
*the standard output.
@@ -67,21 +69,20 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
if (!found) printf("## Particles[???] id=%lld not found\n", id);
}
-void printgParticle(struct gpart *parts, long long int id, size_t N) {
+void printgParticle(struct gpart *gparts, long long int id, size_t N) {
int found = 0;
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
- if (parts[i].id == -id || (parts[i].id > 0 && parts[i].part->id == id)) {
- printf(
- "## gParticle[%zd]: id=%lld, x=[%.16e,%.16e,%.16e], "
- "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], m=%.3e, t_begin=%d, "
- "t_end=%d\n",
- i, parts[i].part->id, parts[i].x[0], parts[i].x[1], parts[i].x[2],
- parts[i].v[0], parts[i].v[1], parts[i].v[2], parts[i].a[0],
- parts[i].a[1], parts[i].a[2], parts[i].mass, parts[i].ti_begin,
- parts[i].ti_end);
+ if (gparts[i].id == -id) {
+ printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
+ gravity_debug_particle(&gparts[i]);
+ found = 1;
+ break;
+ } else if (gparts[i].id > 0 && gparts[i].part->id == id) {
+ printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
+ gravity_debug_particle(&gparts[i]);
found = 1;
break;
}
diff --git a/src/engine.c b/src/engine.c
index 0739609170eceeeae08f1523330e1d7436688a0b..c34214c05b6fb45991208dd78689b58ba5d9731f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -595,7 +595,8 @@ void engine_exchange_cells(struct engine *e) {
* @return The number of arrived parts copied to parts and xparts.
*/
-int engine_exchange_strays(struct engine *e, int offset, size_t *ind, size_t N) {
+int engine_exchange_strays(struct engine *e, int offset, size_t *ind,
+ size_t N) {
#ifdef WITH_MPI
@@ -793,16 +794,18 @@ void engine_maketasks(struct engine *e) {
}
}
- /* Add the gravity mm tasks. */
- for (int i = 0; i < nr_cells; i++)
- if (cells[i].gcount > 0) {
- scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
- &cells[i], NULL, 0);
- for (int j = i + 1; j < nr_cells; j++)
- if (cells[j].gcount > 0)
- scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
- &cells[i], &cells[j], 0);
- }
+ /* /\* Add the gravity mm tasks. *\/ */
+ /* for (int i = 0; i < nr_cells; i++) */
+ /* if (cells[i].gcount > 0) { */
+ /* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1, 0,
+ */
+ /* &cells[i], NULL, 0); */
+ /* for (int j = i + 1; j < nr_cells; j++) */
+ /* if (cells[j].gcount > 0) */
+ /* scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, -1,
+ * 0, */
+ /* &cells[i], &cells[j], 0); */
+ /* } */
/* Split the tasks. */
scheduler_splittasks(sched);
@@ -816,21 +819,24 @@ void engine_maketasks(struct engine *e) {
error("Failed to allocate cell-task links.");
e->nr_links = 0;
- /* Add the gravity up/down tasks at the top-level cells and push them down. */
- for (int k = 0; k < nr_cells; k++)
- if (cells[k].nodeID == nodeID && cells[k].gcount > 0) {
-
- /* Create tasks at top level. */
- struct task *up =
- scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0, 0,
- &cells[k], NULL, 0);
- struct task *down =
- scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0, 0,
- &cells[k], NULL, 0);
-
- /* Push tasks down the cell hierarchy. */
- engine_addtasks_grav(e, &cells[k], up, down);
- }
+ /* /\* Add the gravity up/down tasks at the top-level cells and push them
+ * down. *\/ */
+ /* for (int k = 0; k < nr_cells; k++) */
+ /* if (cells[k].nodeID == nodeID && cells[k].gcount > 0) { */
+
+ /* /\* Create tasks at top level. *\/ */
+ /* struct task *up = */
+ /* scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0,
+ * 0, */
+ /* &cells[k], NULL, 0); */
+ /* struct task *down = */
+ /* scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0,
+ * 0, */
+ /* &cells[k], NULL, 0); */
+
+ /* /\* Push tasks down the cell hierarchy. *\/ */
+ /* engine_addtasks_grav(e, &cells[k], up, down); */
+ /* } */
/* Count the number of tasks associated with each cell and
store the density tasks in each cell, and make each sort
@@ -878,17 +884,17 @@ void engine_maketasks(struct engine *e) {
}
}
- /* Link gravity multipole tasks to the up/down tasks. */
- if (t->type == task_type_grav_mm ||
- (t->type == task_type_sub && t->subtype == task_subtype_grav)) {
- atomic_inc(&t->ci->nr_tasks);
- scheduler_addunlock(sched, t->ci->grav_up, t);
- scheduler_addunlock(sched, t, t->ci->grav_down);
- if (t->cj != NULL && t->ci->grav_up != t->cj->grav_up) {
- scheduler_addunlock(sched, t->cj->grav_up, t);
- scheduler_addunlock(sched, t, t->cj->grav_down);
- }
- }
+ /* /\* Link gravity multipole tasks to the up/down tasks. *\/ */
+ /* if (t->type == task_type_grav_mm || */
+ /* (t->type == task_type_sub && t->subtype == task_subtype_grav)) { */
+ /* atomic_inc(&t->ci->nr_tasks); */
+ /* scheduler_addunlock(sched, t->ci->grav_up, t); */
+ /* scheduler_addunlock(sched, t, t->ci->grav_down); */
+ /* if (t->cj != NULL && t->ci->grav_up != t->cj->grav_up) { */
+ /* scheduler_addunlock(sched, t->cj->grav_up, t); */
+ /* scheduler_addunlock(sched, t, t->cj->grav_down); */
+ /* } */
+ /* } */
}
/* Append a ghost task to each cell, and add kick tasks to the
@@ -1165,8 +1171,8 @@ void engine_print_task_counts(struct engine *e) {
else
counts[task_type_count] += 1;
#ifdef WITH_MPI
- printf("[%04i] %s engine_print_task_counts: task counts are [ %s=%i", e->nodeID,
- clocks_get_timesincestart(), taskID_names[0], counts[0]);
+ printf("[%04i] %s engine_print_task_counts: task counts are [ %s=%i",
+ e->nodeID, clocks_get_timesincestart(), taskID_names[0], counts[0]);
#else
printf("%s engine_print_task_counts: task counts are [ %s=%i",
clocks_get_timesincestart(), taskID_names[0], counts[0]);
diff --git a/src/gravity.h b/src/gravity.h
index d5361095c744e2df069ccaf99447dcd58e10dcd6..f737a0ab882592cffb974f66ccbb62fa3d16d408 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -24,6 +24,6 @@
/* So far only one model here */
/* Straight-forward import */
#include "./gravity/Default/gravity.h"
-#include "./gravity/Default/runner_iact_grav.h"
+#include "./gravity/Default/gravity_iact.h"
#endif
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..98e0c40a5700b4da70f27fb0955592bb5d2287c3
--- /dev/null
+++ b/src/gravity/Default/gravity_debug.h
@@ -0,0 +1,28 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+__attribute__((always_inline))
+ INLINE static void gravity_debug_particle(struct gpart* p) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
+ "mass=%.3e t_begin=%d, t_end=%d\n",
+ p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
+ p->a[0], p->a[1], p->a[2], p->mass, p->ti_begin, p->ti_end);
+}
diff --git a/src/gravity/Default/runner_iact_grav.h b/src/gravity/Default/gravity_iact.h
similarity index 100%
rename from src/gravity/Default/runner_iact_grav.h
rename to src/gravity/Default/gravity_iact.h
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..d707d69631e65eed8ad21a7fa9601c07d3c71263
--- /dev/null
+++ b/src/gravity/Default/gravity_io.h
@@ -0,0 +1,75 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @brief Read dark matter particles from HDF5.
+ *
+ * @param h_grp The HDF5 group in which to read the arrays.
+ * @param N The number of particles on that MPI rank.
+ * @param N_total The total number of particles (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param gparts The particle array
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_read_particles(
+ hid_t h_grp, int N, long long N_total, long long offset,
+ struct gpart* gparts) {
+
+ /* Read arrays */
+ readArray(h_grp, "Coordinates", DOUBLE, N, 3, gparts, N_total, offset, x,
+ COMPULSORY);
+ readArray(h_grp, "Masses", FLOAT, N, 1, gparts, N_total, offset, mass,
+ COMPULSORY);
+ readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v_full,
+ COMPULSORY);
+ readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, gparts, N_total, offset, id,
+ COMPULSORY);
+}
+
+/**
+ * @brief Writes the different particles to the HDF5 file
+ *
+ * @param h_grp The HDF5 group in which to write the arrays.
+ * @param fileName The name of the file (unsued in MPI mode).
+ * @param xmfFile The XMF file to write to (unused in MPI mode).
+ * @param Ndm The number of DM particles on that MPI rank.
+ * @param Ndm_total The total number of g-particles (only used in MPI mode)
+ * @param mpi_rank The MPI rank of this node (only used in MPI mode)
+ * @param offset The offset of the particles for this MPI rank (only used in MPI
+ *mode)
+ * @param gparts The #gpart array
+ * @param us The unit system to use
+ *
+ */
+__attribute__((always_inline)) INLINE static void darkmatter_write_particles(
+ hid_t h_grp, char* fileName, FILE* xmfFile, int Ndm, long long Ndm_total,
+ int mpi_rank, long long offset, struct gpart* gparts,
+ struct UnitSystem* us) {
+
+ /* Write arrays */
+ writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, Ndm, 3, gparts,
+ Ndm_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+ writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, Ndm, 1, gparts,
+ Ndm_total, mpi_rank, offset, mass, us, UNIT_CONV_MASS);
+ writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, Ndm, 3, gparts,
+ Ndm_total, mpi_rank, offset, v_full, us, UNIT_CONV_SPEED);
+ writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, Ndm, 1, gparts,
+ Ndm_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 17262f9d413b1e801b71b456eb4e6baeb925febc..7ce7b81892582f2a90f7dd07f7f244c0d4ed8afb 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -26,7 +26,7 @@ struct gpart {
double x[3];
/* Particle velocity. */
- float v[3];
+ float v_full[3];
/* Particle acceleration. */
float a[3];
@@ -44,7 +44,7 @@ struct gpart {
union {
/* Particle ID. */
- size_t id;
+ long long id;
/* Pointer to corresponding SPH part. */
struct part* part;
diff --git a/src/gravity_io.h b/src/gravity_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..6276e50473de64abd1014bdb36a63a14e02ca8cf
--- /dev/null
+++ b/src/gravity_io.h
@@ -0,0 +1,26 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GRAVITY_IO_H
+#define SWIFT_GRAVITY_IO_H
+
+#include "./const.h"
+
+#include "./gravity/Default/gravity_io.h"
+
+#endif /* SWIFT_GRAVITY_IO_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index 4dbc0388ac20cd7f354007e2e80aa6fccccdc5ad..60453d0c7995f7af2a3166502a24aa590873a043 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -111,7 +111,7 @@ struct part {
float mass;
/* Particle ID. */
- unsigned long long id;
+ long long id;
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 298cfc5ac1d801f03f7ba104a819fd3de531d55c..05754d07dd70bed071e99c86b95eb17eb2194012 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -104,7 +104,7 @@ struct part {
float div_v;
/* Particle ID. */
- unsigned long long id;
+ long long id;
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index 8875e737c6f3c0ab0b7a390b56e9b946c145fc1b..173397ef2c72ee99f4d10742f3645afd1e706218 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -101,7 +101,7 @@ struct part {
} force;
};
- unsigned long long id; /*!< Particle unique ID. */
+ long long id; /*!< Particle unique ID. */
struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
diff --git a/src/map.c b/src/map.c
index dcaa53465767d414cc54fe05940069ae5ff06d77..da13fbfb4ac00ed58184f7fe818826c82265a1de 100644
--- a/src/map.c
+++ b/src/map.c
@@ -96,14 +96,12 @@ void map_cells_plot(struct cell *c, void *data) {
void map_cellcheck(struct cell *c, void *data) {
- int k, *count = (int *)data;
- struct part *p;
-
+ int *count = (int *)data;
__sync_fetch_and_add(count, c->count);
/* Loop over all parts and check if they are in the cell. */
- for (k = 0; k < c->count; k++) {
- p = &c->parts[k];
+ for (int k = 0; k < c->count; k++) {
+ struct part *p = &c->parts[k];
if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
p->x[2] > c->loc[2] + c->h[2]) {
@@ -115,6 +113,21 @@ void map_cellcheck(struct cell *c, void *data) {
error("particle out of bounds!");
}
}
+
+ /* Loop over all gparts and check if they are in the cell. */
+ for (int k = 0; k < c->gcount; k++) {
+ struct gpart *p = &c->gparts[k];
+ if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
+ p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
+ p->x[2] > c->loc[2] + c->h[2]) {
+ printf(
+ "map_cellcheck: g-particle at [ %.16e %.16e %.16e ] outside of cell "
+ "[ %.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n",
+ p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
+ c->loc[0] + c->h[0], c->loc[1] + c->h[1], c->loc[2] + c->h[2]);
+ error("g-particle out of bounds!");
+ }
+ }
}
/**
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 70650d48020cd75eedfad121804dd91923660b08..cffa99a0fd75566ec3e850076d15e104504eeb40 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -357,7 +357,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
*
*/
void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
- int* N, int* periodic, int mpi_rank, int mpi_size,
+ size_t* 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] = {
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 5ee2714a18dbabdd5fc8c07c0b444e04a87ad4f8..a0589944ec845c712abde1e64e305980748db0e7 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -32,7 +32,7 @@
#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
- int* N, int* periodic, int mpi_rank, int mpi_size,
+ size_t* N, int* periodic, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info);
void write_output_parallel(struct engine* e, struct UnitSystem* us,
diff --git a/src/part.h b/src/part.h
index 168e80b68bb211d1d773f1f80f1fa9cc757edfcb..865403e8c2c157dc5a8ff7a32bc41be676d7919b 100644
--- a/src/part.h
+++ b/src/part.h
@@ -35,6 +35,7 @@
/* Some constants. */
#define part_align 64
+#define gpart_align 32
#define xpart_align 32
/* Import the right particle definition */
diff --git a/src/scheduler.c b/src/scheduler.c
index 6c4c7c75a67dc775d04fbddb5e2b479af716963d..722e344b5a86b5fbdc42c7038fd3cb00e44b2ee8 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -534,7 +534,7 @@ void scheduler_splittasks(struct scheduler *s) {
else if (ci->split) {
/* Make a single sub-task? */
- if (scheduler_dosub && ci->count < space_subsize / ci->count) {
+ if (scheduler_dosub && ci->gcount < space_subsize / ci->gcount) {
t->type = task_type_sub;
t->subtype = task_subtype_grav;
@@ -582,7 +582,7 @@ void scheduler_splittasks(struct scheduler *s) {
else {
/* Make a sub-task? */
- if (scheduler_dosub && ci->count < space_subsize / cj->count) {
+ if (scheduler_dosub && ci->gcount < space_subsize / cj->gcount) {
t->type = task_type_sub;
t->subtype = task_subtype_grav;
diff --git a/src/serial_io.c b/src/serial_io.c
index b56d8668b6f4712890120ccbca6e322565ad1775..8e63db5cfad3a3b50fc7e350bbac6ce09708230a 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -410,9 +410,9 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* Calls #error() if an error occurs.
*
*/
-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) {
+void read_ic_serial(char* fileName, double dim[3], struct part** parts,
+ size_t* 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) */
diff --git a/src/serial_io.h b/src/serial_io.h
index e5ecca9c8cbafbbaf2e555c5c216b494e25cc922..95f09f5977a97a359e978db7a1b71b02030d6a14 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -31,9 +31,9 @@
#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
-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);
+void read_ic_serial(char* fileName, double dim[3], struct part** parts,
+ size_t* N, 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 f3225ecf487d223c3bea560cbae626a03da36928..59686a68b5d9e5ea41267ba7b3aad9391862fae4 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -39,6 +39,9 @@
#include "common_io.h"
#include "error.h"
+#define FILENAME_BUFFER_SIZE 150
+#define PARTICLE_GROUP_BUFFER_SIZE 20
+
/*-----------------------------------------------------------------------------
* Routines reading an IC file
*-----------------------------------------------------------------------------*/
@@ -166,7 +169,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
char* temp_c = 0;
hsize_t shape[2];
hsize_t chunk_shape[2];
- char buffer[150];
+ char buffer[FILENAME_BUFFER_SIZE];
/* message("Writing '%s' array...", name); */
@@ -201,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) {
@@ -302,14 +305,18 @@ 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)
*
* @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 dim (output) The dimension of the volume.
+ * @param parts (output) Array of Gas particles.
+ * @param gparts (output) Array of DM particles.
+ * @param Ngas (output) number of Gas particles read.
+ * @param Ngparts (output) The number of DM particles read.
* @param periodic (output) 1 if the volume is periodic, 0 if not.
*
* Opens the HDF5 file fileName and reads the particles contained
@@ -322,13 +329,15 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* Calls #error() if an error occurs.
*
*/
-void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
+void read_ic_single(char* fileName, double dim[3], struct part** parts,
+ struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
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 6 particle types. We only keep the type 0 & 1 for now...*/
+ int numParticles[NUM_PARTICLE_TYPES] = {0};
+ size_t Ndm;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -357,35 +366,82 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
readAttribute(h_grp, "BoxSize", DOUBLE, boxSize);
readAttribute(h_grp, "NumPart_Total", UINT, numParticles);
- *N = numParticles[0];
+ *Ngas = numParticles[0];
+ Ndm = numParticles[1];
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]); */
+ /* *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
/* Close header */
H5Gclose(h_grp);
- /* Allocate memory to store particles */
- 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));
+ /* Total number of particles */
+ *Ngparts = *Ngas + Ndm;
+
+ /* Allocate memory to store SPH particles */
+ 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 */
+ 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.)); */
/* 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.\n");
+ 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++) {
+
+ /* 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 = 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 */
- hydro_read_particles(h_grp, *N, *N, 0, *parts);
+ /* message("Group %s found - reading...", partTypeGroupName); */
- /* Close particle group */
- H5Gclose(h_grp);
+ /* Read particle fields into the particle structure */
+ switch (ptype) {
+
+ case GAS:
+ hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
+ break;
+
+ case DM:
+ darkmatter_read_particles(h_grp, Ndm, Ndm, 0, *gparts);
+ break;
+
+ default:
+ error("Particle Type %d not yet supported. Aborting", ptype);
+ }
+
+ /* Close particle group */
+ H5Gclose(h_grp);
+ }
+
+ /* 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..."); */
@@ -409,28 +465,35 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
*/
void write_output_single(struct engine* e, struct UnitSystem* us) {
- hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
- int N = e->s->nr_parts;
+ hid_t h_file = 0, h_grp = 0;
+ const size_t Ngas = e->s->nr_parts;
+ const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
- int numParticles[6] = {N, 0};
- int numParticlesHighWord[6] = {0};
int numFiles = 1;
struct part* parts = e->s->parts;
- FILE* xmfFile = 0;
+ struct gpart* gparts = e->s->gparts;
+ struct gpart* dmparts = NULL;
static int outputCount = 0;
+ /* 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};
+
/* 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);
/* First time, we need to create the XMF file */
if (outputCount == 0) createXMFfile();
/* Prepare the XMF file for the new entry */
+ FILE* xmfFile = 0;
xmfFile = prepareXMFfile();
/* Write the part corresponding to this specific output */
- writeXMFheader(xmfFile, N, fileName, e->time);
+ writeXMFheader(xmfFile, Ngas, fileName, e->time);
/* Open file */
/* message("Opening file '%s'.", fileName); */
@@ -458,17 +521,20 @@ 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, 6);
+ 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 */
- writeAttribute(h_grp, "NumPart_Total", UINT, numParticles, 6);
+ writeAttribute(h_grp, "NumPart_Total", UINT, numParticles,
+ NUM_PARTICLE_TYPES);
writeAttribute(h_grp, "NumPart_Total_HighWord", UINT, numParticlesHighWord,
- 6);
- 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);
+ NUM_PARTICLE_TYPES);
+ double MassTable[NUM_PARTICLE_TYPES] = {0., 0., 0., 0., 0., 0.};
+ writeAttribute(h_grp, "MassTable", DOUBLE, MassTable, NUM_PARTICLE_TYPES);
+ writeAttribute(h_grp, "Flag_Entropy_ICs", UINT, numParticlesHighWord,
+ NUM_PARTICLE_TYPES);
writeAttribute(h_grp, "NumFilesPerSnapshot", INT, &numFiles, 1);
/* Close header */
@@ -478,25 +544,65 @@ void write_output_single(struct engine* e, struct UnitSystem* us) {
writeCodeDescription(h_file);
/* Print the SPH parameters */
- h_grpsph = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grpsph < 0) error("Error while creating SPH group");
- writeSPHflavour(h_grpsph);
- H5Gclose(h_grpsph);
+ h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grp);
+ H5Gclose(h_grp);
/* 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++) {
- /* Write particle fields from the particle structure */
- hydro_write_particles(h_grp, fileName, xmfFile, N, N, 0, 0, parts, us);
+ /* Don't do anything if no particle of this kind */
+ if (numParticles[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 = H5Gcreate(h_file, partTypeGroupName, H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) {
+ error("Error while creating particle group.\n");
+ }
+
+ /* message("Writing particle arrays..."); */
+
+ /* Write particle fields from the particle structure */
+ switch (ptype) {
+
+ case GAS:
+ hydro_write_particles(h_grp, fileName, xmfFile, Ngas, Ngas, 0, 0, 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, Ndm, Ndm, 0, 0,
+ dmparts, us);
+
+ /* Free temporary array */
+ free(dmparts);
+ break;
+
+ default:
+ error("Particle Type %d not yet supported. Aborting", ptype);
+ }
+
+ /* Close particle group */
+ H5Gclose(h_grp);
+ }
/* Write LXMF file descriptor */
writeXMFfooter(xmfFile);
diff --git a/src/single_io.h b/src/single_io.h
index f6689901106a2cd5d85a873d4047e1d21edd3547..c5250280e82e1801b2a4a6136d404d09093dd0ec 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -26,7 +26,8 @@
#if defined(HAVE_HDF5) && !defined(WITH_MPI)
-void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
+void read_ic_single(char* fileName, double dim[3], struct part** parts,
+ struct gpart** gparts, size_t* Ngas, size_t* Ndm,
int* periodic);
void write_output_single(struct engine* e, struct UnitSystem* us);
diff --git a/src/space.c b/src/space.c
index 5938101d94b17097277d7fc06a1cd4d3c1181727..fa2541c2872c38f5000e78ef2802d9a4f719f9fc 100644
--- a/src/space.c
+++ b/src/space.c
@@ -308,11 +308,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
void space_rebuild(struct space *s, double cell_max, int verbose) {
- int j, k, cdim[3], nr_parts = s->nr_parts, nr_gparts = s->nr_gparts;
- struct cell *restrict c, *restrict cells;
- struct part *restrict p;
- size_t *ind;
- double ih[3], dim[3];
ticks tic = getticks();
/* Be verbose about this. */
@@ -320,13 +315,13 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Re-grid if necessary, or just re-set the cell data. */
space_regrid(s, cell_max, verbose);
- cells = s->cells;
- /* Run through the particles and get their cell index. */
- // tic = getticks();
- const size_t ind_size = s->size_parts;
- if ((ind = (size_t *)malloc(sizeof(size_t) * ind_size)) == NULL)
- error("Failed to allocate temporary particle indices.");
+ int nr_parts = s->nr_parts;
+ int nr_gparts = s->nr_gparts;
+ struct cell *restrict cells = s->cells;
+
+ double ih[3], dim[3];
+ int cdim[3];
ih[0] = s->ih[0];
ih[1] = s->ih[1];
ih[2] = s->ih[2];
@@ -336,9 +331,16 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
cdim[0] = s->cdim[0];
cdim[1] = s->cdim[1];
cdim[2] = s->cdim[2];
- for (k = 0; k < nr_parts; k++) {
- p = &s->parts[k];
- for (j = 0; j < 3; j++)
+
+ /* Run through the particles and get their cell index. */
+ // tic = getticks();
+ const size_t ind_size = s->size_parts;
+ size_t *ind;
+ if ((ind = (size_t *)malloc(sizeof(size_t) * ind_size)) == NULL)
+ error("Failed to allocate temporary particle indices.");
+ for (int k = 0; k < nr_parts; k++) {
+ struct part *restrict p = &s->parts[k];
+ for (int j = 0; j < 3; j++)
if (p->x[j] < 0.0)
p->x[j] += dim[j];
else if (p->x[j] >= dim[j])
@@ -352,8 +354,8 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
#ifdef WITH_MPI
/* Move non-local parts to the end of the list. */
- int nodeID = s->e->nodeID;
- for (k = 0; k < nr_parts; k++)
+ const int nodeID = s->e->nodeID;
+ for (int k = 0; k < nr_parts; k++)
if (cells[ind[k]].nodeID != nodeID) {
cells[ind[k]].count -= 1;
nr_parts -= 1;
@@ -385,8 +387,8 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
}
/* Assign each particle to its cell. */
- for (k = nr_parts; k < s->nr_parts; k++) {
- p = &s->parts[k];
+ for (int k = nr_parts; k < s->nr_parts; k++) {
+ struct part *p = &s->parts[k];
ind[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
cells[ind[k]].count += 1;
@@ -401,7 +403,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1, verbose);
/* Re-link the gparts. */
- for (k = 0; k < nr_parts; k++)
+ for (int k = 0; k < nr_parts; k++)
if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
/* Verify space_sort_struct. */
@@ -419,41 +421,88 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Run through the gravity particles and get their cell index. */
// tic = getticks();
- if ((ind = (size_t *)malloc(sizeof(size_t) * s->size_gparts)) == NULL)
- error("Failed to allocate temporary particle indices.");
- for (k = 0; k < nr_gparts; k++) {
+ const size_t gind_size = s->size_gparts;
+ size_t *gind;
+ if ((gind = (size_t *)malloc(sizeof(size_t) * gind_size)) == NULL)
+ error("Failed to allocate temporary g-particle indices.");
+ for (int k = 0; k < nr_gparts; k++) {
struct gpart *gp = &s->gparts[k];
- for (j = 0; j < 3; j++)
+ for (int j = 0; j < 3; j++)
if (gp->x[j] < 0.0)
gp->x[j] += dim[j];
else if (gp->x[j] >= dim[j])
gp->x[j] -= dim[j];
- ind[k] =
+ gind[k] =
cell_getid(cdim, gp->x[0] * ih[0], gp->x[1] * ih[1], gp->x[2] * ih[2]);
- cells[ind[k]].gcount++;
+ cells[gind[k]].gcount++;
}
- // message( "getting particle indices took %.3f %s." ,
- // clocks_from_ticks(getticks() - tic), clocks_getunit());
+// message( "getting particle indices took %.3f %s." ,
+// clocks_from_ticks(getticks() - tic), clocks_getunit());
+
+#ifdef WITH_MPI
+
+ /* Move non-local gparts to the end of the list. */
+ for (int k = 0; k < nr_gparts; k++)
+ if (cells[ind[k]].nodeID != nodeID) {
+ cells[ind[k]].gcount -= 1;
+ nr_gparts -= 1;
+ struct gpart tp = s->gparts[k];
+ s->gparts[k] = s->gparts[nr_gparts];
+ s->gparts[nr_gparts] = tp;
+ int t = ind[k];
+ ind[k] = ind[nr_gparts];
+ ind[nr_gparts] = t;
+ }
+
+ /* Exchange the strays, note that this potentially re-allocates
+ the parts arrays. */
+ // s->nr_gparts =
+ // nr_gparts + engine_exchange_strays(s->e, nr_gparts, &ind[nr_gparts],
+ // s->nr_gparts - nr_gparts);
+ if (nr_gparts > 0)
+ error("Need to implement the exchange of strays for the gparts");
+
+ /* Re-allocate the index array if needed.. */
+ if (s->nr_gparts > gind_size) {
+ size_t *gind_new;
+ if ((gind_new = (size_t *)malloc(sizeof(size_t) * s->nr_gparts)) == NULL)
+ error("Failed to allocate temporary g-particle indices.");
+ memcpy(gind_new, gind, sizeof(size_t) * nr_gparts);
+ free(gind);
+ gind = gind_new;
+ }
+
+ /* Assign each particle to its cell. */
+ for (int k = nr_gparts; k < s->nr_gparts; k++) {
+ struct gpart *p = &s->gparts[k];
+ gind[k] =
+ cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
+ cells[gind[k]].count += 1;
+ /* if ( cells[ ind[k] ].nodeID != nodeID )
+ error( "Received part that does not belong to me (nodeID=%i)." , cells[
+ ind[k] ].nodeID ); */
+ }
+ nr_gparts = s->nr_gparts;
- /* TODO: Here we should exchange the gparts as well! */
+#endif
/* Sort the parts according to their cells. */
- space_gparts_sort(s->gparts, ind, nr_gparts, 0, s->nr_cells - 1);
+ space_gparts_sort(s->gparts, gind, nr_gparts, 0, s->nr_cells - 1);
/* Re-link the parts. */
- for (k = 0; k < nr_gparts; k++)
+ for (int k = 0; k < nr_gparts; k++)
if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
/* We no longer need the indices as of here. */
- free(ind);
+ free(gind);
/* Hook the cells up to the parts. */
// tic = getticks();
struct part *finger = s->parts;
struct xpart *xfinger = s->xparts;
struct gpart *gfinger = s->gparts;
- for (k = 0; k < s->nr_cells; k++) {
- c = &cells[k];
+ for (int k = 0; k < s->nr_cells; k++) {
+ struct cell *restrict c = &cells[k];
c->parts = finger;
c->xparts = xfinger;
c->gparts = gfinger;
@@ -1003,12 +1052,15 @@ void space_map_cells_pre(struct space *s, int full,
void space_do_split(struct space *s, struct cell *c) {
- int k, count = c->count, gcount = c->gcount, maxdepth = 0;
- float h, h_max = 0.0f;
- int ti_end_min = max_nr_timesteps, ti_end_max = 0, ti_end;
+ const int count = c->count;
+ const int gcount = c->gcount;
+ int maxdepth = 0;
+ float h_max = 0.0f;
+ int ti_end_min = max_nr_timesteps, ti_end_max = 0;
struct cell *temp;
- struct part *p, *parts = c->parts;
- struct xpart *xp, *xparts = c->xparts;
+ struct part *parts = c->parts;
+ struct gpart *gparts = c->gparts;
+ struct xpart *xparts = c->xparts;
/* Check the depth. */
if (c->depth > s->maxdepth) s->maxdepth = c->depth;
@@ -1020,7 +1072,7 @@ void space_do_split(struct space *s, struct cell *c) {
c->split = 1;
/* Create the cell's progeny. */
- for (k = 0; k < 8; k++) {
+ for (int k = 0; k < 8; k++) {
temp = space_getcell(s);
temp->count = 0;
temp->gcount = 0;
@@ -1047,7 +1099,7 @@ void space_do_split(struct space *s, struct cell *c) {
cell_split(c);
/* Remove any progeny with zero parts. */
- for (k = 0; k < 8; k++)
+ for (int k = 0; k < 8; k++)
if (c->progeny[k]->count == 0 && c->progeny[k]->gcount == 0) {
space_recycle(s, c->progeny[k]);
c->progeny[k] = NULL;
@@ -1077,26 +1129,36 @@ void space_do_split(struct space *s, struct cell *c) {
c->maxdepth = c->depth;
/* Get dt_min/dt_max. */
-
- for (k = 0; k < count; k++) {
- p = &parts[k];
- xp = &xparts[k];
+ for (int k = 0; k < count; k++) {
+ struct part *p = &parts[k];
+ struct xpart *xp = &xparts[k];
+ const float h = p->h;
+ const int ti_end = p->ti_end;
xp->x_old[0] = p->x[0];
xp->x_old[1] = p->x[1];
xp->x_old[2] = p->x[2];
- h = p->h;
- ti_end = p->ti_end;
if (h > h_max) h_max = h;
if (ti_end < ti_end_min) ti_end_min = ti_end;
if (ti_end > ti_end_max) ti_end_max = ti_end;
}
+ for (int k = 0; k < gcount; k++) {
+ struct gpart *p = &gparts[k];
+ const int ti_end = p->ti_end;
+ if (ti_end < ti_end_min) ti_end_min = ti_end;
+ if (ti_end > ti_end_max) ti_end_max = ti_end;
+ }
c->h_max = h_max;
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
}
/* Set ownership according to the start of the parts array. */
- c->owner = ((c->parts - s->parts) % s->nr_parts) * s->nr_queues / s->nr_parts;
+ if (s->nr_parts > 0)
+ c->owner =
+ ((c->parts - s->parts) % s->nr_parts) * s->nr_queues / s->nr_parts;
+ else
+ c->owner =
+ ((c->gparts - s->gparts) % s->nr_gparts) * s->nr_queues / s->nr_gparts;
}
/**
@@ -1176,8 +1238,10 @@ struct cell *space_getcell(struct space *s) {
*
* @param s The #space to initialize.
* @param dim Spatial dimensions of the domain.
- * @param parts Pointer to an array of #part.
- * @param N The number of parts in the space.
+ * @param parts Array of Gas particles.
+ * @param gparts Array of Gravity particles.
+ * @param Ngas The number of Gas particles in the space.
+ * @param Ngpart The number of Gravity particles in the space.
* @param periodic flag whether the domain is periodic or not.
* @param h_max The maximal interaction radius.
* @param verbose Print messages to stdout or not
@@ -1188,62 +1252,59 @@ struct cell *space_getcell(struct space *s) {
* recursively.
*/
-void space_init(struct space *s, double dim[3], struct part *parts, size_t N,
- int periodic, double h_max, int verbose) {
+void space_init(struct space *s, double dim[3], struct part *parts,
+ struct gpart *gparts, size_t Ngas, size_t Ngpart, int periodic,
+ double h_max, int verbose) {
/* Store everything in the space. */
s->dim[0] = dim[0];
s->dim[1] = dim[1];
s->dim[2] = dim[2];
s->periodic = periodic;
- s->nr_parts = N;
- s->size_parts = N;
+ s->nr_parts = Ngas;
+ s->size_parts = Ngas;
s->parts = parts;
+ s->nr_gparts = Ngpart;
+ s->size_gparts = Ngpart;
+ s->gparts = gparts;
s->cell_min = h_max;
s->nr_queues = 1;
s->size_parts_foreign = 0;
- /* Check that all the particle positions are reasonable, wrap if periodic. */
+ /* Check that all the gas particle positions are reasonable, wrap if periodic.
+ */
if (periodic) {
- for (int k = 0; k < N; k++)
+ for (int k = 0; k < Ngas; k++)
for (int j = 0; j < 3; j++) {
while (parts[k].x[j] < 0) parts[k].x[j] += dim[j];
while (parts[k].x[j] >= dim[j]) parts[k].x[j] -= dim[j];
}
} else {
- for (int k = 0; k < N; k++)
+ for (int k = 0; k < Ngas; k++)
for (int j = 0; j < 3; j++)
if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
error("Not all particles are within the specified domain.");
}
- /* Allocate the xtra parts array. */
- if (posix_memalign((void *)&s->xparts, part_align,
- N * sizeof(struct xpart)) != 0)
- error("Failed to allocate xparts.");
- bzero(s->xparts, N * sizeof(struct xpart));
-
- /* For now, clone the parts to make gparts. */
- if (posix_memalign((void *)&s->gparts, part_align,
- N * sizeof(struct gpart)) != 0)
- error("Failed to allocate gparts.");
- bzero(s->gparts, N * sizeof(struct gpart));
- /* for ( int k = 0 ; k < N ; k++ ) {
- s->gparts[k].x[0] = s->parts[k].x[0];
- s->gparts[k].x[1] = s->parts[k].x[1];
- s->gparts[k].x[2] = s->parts[k].x[2];
- s->gparts[k].v[0] = s->parts[k].v[0];
- s->gparts[k].v[1] = s->parts[k].v[1];
- s->gparts[k].v[2] = s->parts[k].v[2];
- s->gparts[k].mass = s->parts[k].mass;
- s->gparts[k].dt = s->parts[k].dt;
- s->gparts[k].id = s->parts[k].id;
- s->gparts[k].part = &s->parts[k];
- s->parts[k].gpart = &s->gparts[k];
+ /* Same for the gparts */
+ if (periodic) {
+ for (int k = 0; k < Ngpart; k++)
+ for (int j = 0; j < 3; j++) {
+ while (gparts[k].x[j] < 0) gparts[k].x[j] += dim[j];
+ while (gparts[k].x[j] >= dim[j]) gparts[k].x[j] -= dim[j];
}
- s->nr_gparts = s->nr_parts; */
- s->nr_gparts = 0;
- s->size_gparts = s->size_parts;
+ } else {
+ for (int k = 0; k < Ngpart; k++)
+ for (int j = 0; j < 3; j++)
+ if (gparts[k].x[j] < 0 || gparts[k].x[j] >= dim[j])
+ error("Not all particles are within the specified domain.");
+ }
+
+ /* Allocate the extra parts array. */
+ if (posix_memalign((void *)&s->xparts, xpart_align,
+ Ngas * sizeof(struct xpart)) != 0)
+ error("Failed to allocate xparts.");
+ bzero(s->xparts, Ngas * sizeof(struct xpart));
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
diff --git a/src/space.h b/src/space.h
index 87b618b141b7871aeb5ff76eae5bdc24f063bdef..91485ff7e2ebe9da8ab927748589ae9f71320803 100644
--- a/src/space.h
+++ b/src/space.h
@@ -127,12 +127,14 @@ extern struct parallel_sort space_sort_struct;
/* function prototypes. */
void space_parts_sort(struct space *s, size_t *ind, size_t N, int min, int max,
int verbose);
-void space_gparts_sort(struct gpart *gparts, size_t *ind, size_t N, int min, int max);
+void space_gparts_sort(struct gpart *gparts, size_t *ind, size_t N, int min,
+ int max);
struct cell *space_getcell(struct space *s);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
-void space_init(struct space *s, double dim[3], struct part *parts, size_t N,
- int periodic, double h_max, int verbose);
+void space_init(struct space *s, double dim[3], struct part *parts,
+ struct gpart *gparts, size_t N, size_t Ngpart, int periodic,
+ double h_max, int verbose);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
void space_map_parts(struct space *s,
diff --git a/tests/testReading.c b/tests/testReading.c
index 17b33d88f6fd702e0bad6280f4adbe7ebc50ef2a..d2a2a766171a85ace486914f0f39a987d9d8c3d3 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -22,10 +22,12 @@
int main() {
- int N = -1, periodic = -1;
+ int Ngas = -1, Ngpart = -1;
+ int periodic = -1;
int i, j, k, n;
double dim[3];
struct part *parts = NULL;
+ struct gpart *gparts = NULL;
/* Properties of the ICs */
const double boxSize = 1.;
@@ -33,25 +35,26 @@ int main() {
const double rho = 2.;
/* Read data */
- read_ic_single("input.hdf5", dim, &parts, &N, &periodic);
+ read_ic_single("input.hdf5", dim, &parts, &gparts, &Ngas, &Ngpart, &periodic);
/* Check global properties read are correct */
assert(dim[0] == boxSize);
assert(dim[1] == boxSize);
assert(dim[2] == boxSize);
- assert(N == L * L * L);
+ assert(Ngas == L * L * L);
+ assert(Ngpart == L * L * L);
assert(periodic == 1);
/* Check particles */
- for (n = 0; n < N; ++n) {
+ for (n = 0; n < Ngas; ++n) {
/* Check that indices are in a reasonable range */
unsigned long long index = parts[n].id;
- assert(index < N);
+ assert(index < Ngas);
/* Check masses */
float mass = parts[n].mass;
- float correct_mass = boxSize * boxSize * boxSize * rho / N;
+ float correct_mass = boxSize * boxSize * boxSize * rho / Ngas;
assert(mass == correct_mass);
/* Check smoothing length */