From ee225bf7436587b58fa4bdba6cda2b7ed85c6508 Mon Sep 17 00:00:00 2001
From: John Regan <john.a.regan@durham.ac.uk>
Date: Fri, 19 Feb 2016 11:10:44 +0000
Subject: [PATCH] First additions to the GravityParticles branch. IO working,
everything else broken
---
examples/main.c | 36 +++++---
src/common_io.h | 8 ++
src/engine.c | 6 +-
src/gravity/Default/gravity_io.h | 71 +++++++++++++++
src/gravity_io.h | 26 ++++++
src/single_io.c | 70 +++++++++++----
src/single_io.h | 4 +-
src/space.c | 147 ++++++++++++++++---------------
src/space.h | 4 +-
9 files changed, 267 insertions(+), 105 deletions(-)
create mode 100644 src/gravity/Default/gravity_io.h
create mode 100644 src/gravity_io.h
diff --git a/examples/main.c b/examples/main.c
index 32264f185f..07bf70e2e6 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -63,7 +63,7 @@
int main(int argc, char *argv[]) {
- int c, icount, j, k, N, periodic = 1;
+ int c, icount, j, k, Ngas, Ndm, periodic = 1;
long long N_total = -1;
int nr_threads = 1, nr_queues = -1;
int dump_tasks = 0;
@@ -72,6 +72,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;
@@ -264,14 +265,14 @@ int main(int argc, char *argv[]) {
tic = getticks();
#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, &gparts, &Ngas, &Ndm, &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, &gparts, &Ngas, &Ndm, &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, &Ndm, &periodic);
#endif
if (myrank == 0)
@@ -280,31 +281,41 @@ int main(int argc, char *argv[]) {
fflush(stdout);
#if defined(WITH_MPI)
- long long N_long = N;
+ long long tmp;
+ long long N_long = Ngas;
+ MPI_Reduce(&N_long, &tmp, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+ long long N_long = Ndm;
MPI_Reduce(&N_long, &N_total, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
+ N_total += tmp;
#else
- N_total = N;
+ N_total = Ngas + Ndm;
#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 (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 (k = 0; k < Ngas; k++) {
parts[k].x[0] += shift[0];
parts[k].x[1] += shift[1];
parts[k].x[2] += shift[2];
}
+ for (k = 0; k < Ndm; 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;
/* Initialize the space with this data. */
tic = getticks();
- space_init(&s, dim, parts, N, periodic, h_max, myrank == 0);
+ space_init(&s, dim, parts, gparts, Ngas, Ndm, periodic, h_max, myrank == 0);
if (myrank == 0)
message("space_init took %.3f ms.",
((double)(getticks() - tic)) / CPU_TPS * 1000);
@@ -317,7 +328,8 @@ int main(int argc, char *argv[]) {
message("space %s periodic.", s.periodic ? "is" : "isn't");
message("highest-level cell dimensions are [ %i %i %i ].", s.cdim[0],
s.cdim[1], s.cdim[2]);
- message("%i parts in %i cells.", s.nr_parts, s.tot_cells);
+ message("%i gas parts in %i cells.", s.nr_parts, s.tot_cells);
+ message("%i dm parts 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);
}
@@ -389,7 +401,7 @@ int main(int argc, char *argv[]) {
/* Initialise the particles */
engine_init_particles(&e);
-
+ exit(-99);
/* Legend */
if (myrank == 0)
printf(
diff --git a/src/common_io.h b/src/common_io.h
index c221ad3aaf..1447d291ea 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -55,6 +55,14 @@ enum DATA_IMPORTANCE {
OPTIONAL = 0
};
+#define NUMBER_PARTICLE_TYPES 6
+#define PARTICLE_TYPE_STRINGLEN 12
+
+/* Particle Types can contain up to NUMBER_PARTICLE_TYPES types. */
+enum PARTICLE_TYPE {
+ GAS = 0,
+ DARKMATTER
+};
hid_t hdf5Type(enum DATA_TYPE type);
size_t sizeOfType(enum DATA_TYPE type);
diff --git a/src/engine.c b/src/engine.c
index 2cca53610c..413cfb4a2e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1549,7 +1549,7 @@ void engine_prepare(struct engine *e) {
// tic)/CPU_TPS*1000 );
#endif
e->tic_step = getticks();
-
+
/* Did this not go through? */
if (rebuild) {
// tic = getticks();
@@ -1557,7 +1557,7 @@ void engine_prepare(struct engine *e) {
// message( "engine_rebuild took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
}
-
+
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
// tic = getticks();
@@ -1718,7 +1718,7 @@ void engine_init_particles(struct engine *e) {
struct space *s = e->s;
message("Initialising particles");
-
+
engine_prepare(e);
/* Make sure all particles are ready to go */
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
new file mode 100644
index 0000000000..796d24a93d
--- /dev/null
+++ b/src/gravity/Default/gravity_io.h
@@ -0,0 +1,71 @@
+/*******************************************************************************
+ * 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 Reads the different particles to the HDF5 file
+ *
+ * @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 parts 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) {
+ message("Reading Dark Matter particles\n");
+ /* Read arrays */
+ readArray(h_grp, "Coordinates", DOUBLE, N, 3, gparts, N_total, offset, x,
+ COMPULSORY);
+ readArray(h_grp, "Velocities", FLOAT, N, 3, gparts, N_total, offset, v,
+ 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 N The number of particles on that MPI rank.
+ * @param N_total The total number of 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 parts The particle 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 N, long long N_total,
+ int mpi_rank, long long offset, struct gpart* gparts, struct UnitSystem* us) {
+
+ /* Write arrays */
+ writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, gparts,
+ N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
+ writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, gparts,
+ N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
+ writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, gparts,
+ N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
+}
+
diff --git a/src/gravity_io.h b/src/gravity_io.h
new file mode 100644
index 0000000000..6276e50473
--- /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/single_io.c b/src/single_io.c
index b0e350ffe2..2a03554a4e 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -276,6 +276,8 @@ 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)
@@ -296,13 +298,20 @@ 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,
- int* periodic) {
+void read_ic_single(char* fileName, double dim[3], struct part** parts, struct gpart** gparts,
+ int *Ngas, int *Ndm, 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*/
+ int numParticles[NUMBER_PARTICLE_TYPES] = {0};
+ int ptype = 0;
+ char partType[PARTICLE_TYPE_STRINGLEN];
+ /* HDF5 parameters */
+ herr_t group_exists = -1;
+ hid_t lapl_id = H5Pcreate(H5P_LINK_ACCESS);
+
+
+ /* GADGET has 6 particle types. We only keep the type 0 & 1 for now...*/
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
@@ -331,7 +340,8 @@ 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];
@@ -343,23 +353,53 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
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));
+ 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 gravity particles */
+ if (posix_memalign((void*)gparts, part_align, *Ndm * sizeof(struct gpart)) != 0)
+ error("Error while allocating memory for gravity particles");
+ bzero(*gparts, *Ndm * 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 = H5Gopen1(h_file, "/PartType0");
- if (h_grp < 0) error("Error while opening particle group.\n");
-
- /* Read particle fields into the particle structure */
- hydro_read_particles(h_grp, *N, *N, 0, *parts);
- /* Close particle group */
- H5Gclose(h_grp);
+ message("BoxSize = %lf\n", dim[0]);
+ message("NumPart_Total = %d\n", *Ndm + *Ngas);
+ fflush(stdout);
+ /* Loop over all particle types */
+ for(ptype = 0; ptype < NUMBER_PARTICLE_TYPES; ptype++) {
+ snprintf(partType, PARTICLE_TYPE_STRINGLEN, "/PartType%d", ptype);
+ group_exists = H5Lexists(h_file, partType,lapl_id);
+
+ if(group_exists) {
+ h_grp = H5Gopen2(h_file, partType, H5P_DEFAULT);
+ if (h_grp < 0) {
+ error("Error while opening particle group.\n");
+ }
+ else {
+ message("Group %s found - reading.", partType);
+ /* Read particle fields into the particle structure */
+ if(GAS == ptype)
+ hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
+ else if(DARKMATTER == ptype)
+ darkmatter_read_particles(h_grp, *Ndm, *Ndm, 0, *gparts);
+ else
+ error("Particle Type %d not yet supported. Aborting", ptype);
+
+
+ /* Close particle group */
+ H5Gclose(h_grp);
+ }
+ }
+ else {
+ //message("Particle Type %d does not exist - skipping", part);
+ continue;
+ }
+ }
/* message("Done Reading particles..."); */
diff --git a/src/single_io.h b/src/single_io.h
index f668990110..8d795d94ee 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -26,8 +26,8 @@
#if defined(HAVE_HDF5) && !defined(WITH_MPI)
-void read_ic_single(char* fileName, double dim[3], struct part** parts, int* N,
- int* periodic);
+void read_ic_single(char* fileName, double dim[3], struct part** parts, struct gpart** gparts,
+ int* Ngas, int* Ndm, int* periodic);
void write_output_single(struct engine* e, struct UnitSystem* us);
diff --git a/src/space.c b/src/space.c
index 60f86bcde8..74d480790f 100644
--- a/src/space.c
+++ b/src/space.c
@@ -164,14 +164,20 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
/* Run through the parts and get the current h_max. */
// tic = getticks();
- if (s->cells != NULL) {
- for (k = 0; k < s->nr_cells; k++) {
- if (s->cells[k].h_max > h_max) h_max = s->cells[k].h_max;
- }
- } else {
- for (k = 0; k < nr_parts; k++) {
- if (s->parts[k].h > h_max) h_max = s->parts[k].h;
+ if(nr_parts) {
+ if (s->cells != NULL) {
+ for (k = 0; k < s->nr_cells; k++) {
+ if (s->cells[k].h_max > h_max) h_max = s->cells[k].h_max;
+ }
+ } else {
+ for (k = 0; k < nr_parts; k++) {
+ if (s->parts[k].h > h_max) h_max = s->parts[k].h;
+ }
+ s->h_max = h_max;
}
+ }
+ else {
+ h_max = s->dim[0]/16.0;
s->h_max = h_max;
}
@@ -308,6 +314,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
struct part *restrict p;
int *ind;
double ih[3], dim[3];
+
// ticks tic;
/* Be verbose about this. */
@@ -316,10 +323,11 @@ 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. */
+
+ /* Run through the SPH particles and get their cell index. */
// tic = getticks();
- const int ind_size = s->size_parts;
+
+ const int ind_size = nr_parts;
if ((ind = (int *)malloc(sizeof(int) * ind_size)) == NULL)
error("Failed to allocate temporary particle indices.");
ih[0] = s->ih[0];
@@ -332,16 +340,17 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
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++)
- if (p->x[j] < 0.0)
- p->x[j] += dim[j];
- else if (p->x[j] >= dim[j])
- p->x[j] -= dim[j];
- ind[k] =
+ p = &s->parts[k];
+ for (j = 0; j < 3; j++)
+ if (p->x[j] < 0.0)
+ p->x[j] += dim[j];
+ else if (p->x[j] >= dim[j])
+ p->x[j] -= dim[j];
+ ind[k] =
cell_getid(cdim, p->x[0] * ih[0], p->x[1] * ih[1], p->x[2] * ih[2]);
- cells[ind[k]].count++;
+ cells[ind[k]].count++;
}
+
// message( "getting particle indices took %.3f ms." , (double)(getticks() -
// tic) / CPU_TPS * 1000 );
@@ -395,26 +404,14 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Sort the parts according to their cells. */
// tic = getticks();
space_parts_sort(s, ind, nr_parts, 0, s->nr_cells - 1);
+
// message( "parts_sort took %.3f ms." , (double)(getticks() - tic) / CPU_TPS
// * 1000 );
- /* Re-link the gparts. */
- for (k = 0; k < nr_parts; k++)
- if (s->parts[k].gpart != NULL) s->parts[k].gpart->part = &s->parts[k];
-
- /* Verify space_sort_struct. */
- /* for ( k = 1 ; k < nr_parts ; k++ ) {
- if ( ind[k-1] > ind[k] ) {
- error( "Sort failed!" );
- }
- else if ( ind[k] != cell_getid( cdim , parts[k].x[0]*ih[0] ,
- parts[k].x[1]*ih[1] , parts[k].x[2]*ih[2] ) )
- error( "Incorrect indices!" );
- } */
-
+
/* We no longer need the indices as of here. */
free(ind);
-
+
/* Run through the gravity particles and get their cell index. */
// tic = getticks();
if ((ind = (int *)malloc(sizeof(int) * s->size_gparts)) == NULL)
@@ -441,13 +438,10 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
// message( "gparts_sort took %.3f ms." , (double)(getticks() - tic) / CPU_TPS
// * 1000 );
- /* Re-link the parts. */
- for (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);
-
+
/* Hook the cells up to the parts. */
// tic = getticks();
struct part *finger = s->parts;
@@ -462,6 +456,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
xfinger = &xfinger[c->count];
gfinger = &gfinger[c->gcount];
}
+
// message( "hooking up cells took %.3f ms." , (double)(getticks() - tic) /
// CPU_TPS * 1000 );
@@ -980,6 +975,7 @@ void space_split(struct space *s, struct cell *c) {
struct cell *temp;
struct part *p, *parts = c->parts;
struct xpart *xp, *xparts = c->xparts;
+ struct gpart *gp, *gparts = c->gparts;
/* Check the depth. */
if (c->depth > s->maxdepth) s->maxdepth = c->depth;
@@ -1013,7 +1009,7 @@ void space_split(struct space *s, struct cell *c) {
temp->parent = c;
c->progeny[k] = temp;
}
-
+
/* Split the cell data. */
cell_split(c);
@@ -1052,6 +1048,7 @@ void space_split(struct space *s, struct cell *c) {
for (k = 0; k < count; k++) {
p = &parts[k];
xp = &xparts[k];
+
xp->x_old[0] = p->x[0];
xp->x_old[1] = p->x[1];
xp->x_old[2] = p->x[2];
@@ -1060,14 +1057,26 @@ void space_split(struct space *s, struct cell *c) {
if (h > h_max) h_max = h;
if (t_end < t_end_min) t_end_min = t_end;
if (t_end > t_end_max) t_end_max = t_end;
+ c->h_max = h_max;
}
- c->h_max = h_max;
+ if (!count) {
+ for (k = 0; k < count; k++) {
+ gp = &gparts[k];
+ t_end = gp->t_end;
+ if (t_end < t_end_min) t_end_min = t_end;
+ if (t_end > t_end_max) t_end_max = t_end;
+ }
+ }
+
c->t_end_min = t_end_min;
c->t_end_max = t_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(count)
+ 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;
}
/**
@@ -1159,62 +1168,58 @@ struct cell *space_getcell(struct space *s) {
* recursively.
*/
-void space_init(struct space *s, double dim[3], struct part *parts, int N,
- int periodic, double h_max, int verbose) {
+void space_init(struct space *s, double dim[3], struct part *parts, struct gpart *gparts,
+ int Ngas, int Ndm, 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;
+ /* SPH */
s->parts = parts;
+ s->nr_parts = Ngas;
+ s->size_parts = Ngas;
s->cell_min = h_max;
s->nr_queues = 1;
+ /* Dark Matter */
+ s->nr_gparts = Ndm;
+ s->size_gparts = Ndm;
+ s->gparts = gparts;
+
s->size_parts_foreign = 0;
+
+ message("s->nr_parts = %d\t s->cell_min = %g\n", s->nr_parts, s->cell_min);
+ message("s->nr_gparts = %d\n", s->nr_gparts);
/* Check that all the 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];
}
+ for (int k = 0; k < Ndm; 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];
+ }
} 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.");
+ for (int k = 0; k < Ndm; 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 xtra parts array. */
if (posix_memalign((void *)&s->xparts, part_align,
- N * sizeof(struct xpart)) != 0)
+ Ngas * 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];
- }
- s->nr_gparts = s->nr_parts; */
- s->nr_gparts = 0;
- s->size_gparts = s->size_parts;
+ 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 91d3d9a4d3..cffacf3ca8 100644
--- a/src/space.h
+++ b/src/space.h
@@ -126,8 +126,8 @@ void gparts_sort(struct gpart *gparts, int *ind, int 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, int N,
- int periodic, double h_max, int verbose);
+void space_init(struct space *s, double dim[3], struct part *parts, struct gpart *gparts,
+ int Ngas, int Ndm, 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,
--
GitLab