Commit ee225bf7 authored by John Regan's avatar John Regan
Browse files

First additions to the GravityParticles branch. IO working, everything else broken

parent 3b2b2be6
......@@ -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(
......
......@@ -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);
......
......@@ -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 */
......
/*******************************************************************************
* 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);
}
/*******************************************************************************
* 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 */
......@@ -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..."); */
......
......@@ -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);
......
......@@ -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++ ) {