Commit c7fa73cf authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Imported the multi-type i/o and added proper gravity debugging routines

parent 3e0e0c69
......@@ -58,8 +58,8 @@
int main(int argc, char *argv[]) {
int c, icount, j, k, N, periodic = 1;
long long N_total = -1;
int c, icount, j, k, Ngas = 0, Ngpart = 0, periodic = 1;
long long N_total[2] = {0};
int nr_threads = 1, nr_queues = -1;
int dump_tasks = 0;
int data[2];
......@@ -67,6 +67,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;
......@@ -361,7 +362,7 @@ int main(int argc, char *argv[]) {
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 +373,34 @@ 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);
N_total[1] -= N_total[0]
#else
N_total = N;
N_total[0] = Ngas;
N_total[1] = Ngpart - Ngas;
#endif
if (myrank == 0) message("Read %lld particles from the ICs", N_total);
if (myrank == 0)
message("Read %lld gas particles and %lld DM particles from the ICs",
N_total[0], N_total[1]);
/* 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 < 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 +410,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),
......@@ -482,9 +494,11 @@ int main(int argc, char *argv[]) {
if (myrank == 0) {
message(
"Running on %lld particles until t=%.3e with %i threads and %i "
"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, time_end, e.nr_threads, e.sched.nr_queues, e.dt_min, e.dt_max);
N_total[0], N_total[1], time_end, e.nr_threads, e.sched.nr_queues,
e.dt_min, e.dt_max);
fflush(stdout);
}
......
......@@ -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];
......
......@@ -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,85 @@ 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, int Ndm) {
/* Let's give all these gparts a negative id */
for (int i = 0; i < Ndm; ++i) {
gparts[i].id = -abs(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, int Ngas,
int Ndm) {
for (int 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, int Ntot, struct gpart* dmparts,
int Ndm) {
int count = 0;
/* Loop over all gparts */
for (int i = 0; i < Ntot; ++i) {
/* And collect the DM ones */
if (gparts[i].id < 0) {
memcpy(&dmparts[count], &gparts[i], sizeof(struct gpart));
count++;
}
}
/* Check that everything is fine */
if (count != Ndm)
error("Collected the wrong number of dm particles (%d vs. %d expected)",
count, Ndm);
}
#endif
......@@ -24,6 +24,7 @@
#include "../config.h"
/* Includes. */
#include "part.h"
#include "units.h"
#if defined(HAVE_HDF5)
......@@ -55,9 +56,26 @@ enum DATA_IMPORTANCE {
OPTIONAL = 0
};
/**
* @brief The different particle types present in a GADGET IC file
*
*/
enum PARTICLE_TYPE {
GAS = 0,
DM = 1,
STAR = 4,
BH = 5
};
hid_t hdf5Type(enum DATA_TYPE type);
size_t sizeOfType(enum DATA_TYPE type);
void collect_dm_gparts(struct gpart* gparts, int Ntot, struct gpart* dmparts,
int Ndm);
void prepare_dm_gparts(struct gpart* gparts, int Ndm);
void duplicate_hydro_gparts(struct part* parts, struct gpart* gparts, int Ngas,
int 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,
......
......@@ -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;
}
......
......@@ -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
......@@ -1165,8 +1166,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]);
......
/*******************************************************************************
* 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);
}
/*******************************************************************************
* 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) {
/* 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);
}
......@@ -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;
......
/*******************************************************************************
* 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 */
......@@ -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;
......
......@@ -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;
......
......@@ -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. */
......
......@@ -35,6 +35,7 @@
/* Some constants. */
#define part_align 64
#define gpart_align 32
#define xpart_align 32
/* Import the right particle definition */
......
......@@ -199,9 +199,6 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
chunk_shape[1] = 0;
}
/* 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,6 +299,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)
......@@ -322,13 +321,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, int* Ngas, int* 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) */
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[6] = {0};
/* GADGET has 6 particle types. We only keep the type 0*/
int Ndm;
/* Open file */
/* message("Opening file '%s' as IC.", fileName); */
......@@ -357,35 +358,81 @@ 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\n", dim[0]);
message("NumPart = [%d, %d] Total = %d\n", *Ngas, Ndm, *Ngparts);
/* Read particle fields into the particle structure */
hydro_read_particles(h_grp, *N, *N, 0, *parts);
/* Loop over all particle types */
for (int ptype = 0; ptype < 6; ptype++) {
/* Close particle group */
H5Gclose(h_grp);
/* Don't do anything if no particle of this kind */
if (numParticles[ptype] == 0) continue;
/* Open the particle group in the file */
char partTypeGroupName[15];
sprintf(partTypeGroupName, "/PartType%d", ptype);
h_grp = H5Gopen(h_file, partTypeGroupName, H5P_DEFAULT);
if (h_grp < 0) {
error("Error while opening particle group %s.", partTypeGroupName);
}
/* message("Group %s found - reading...", partTypeGroupName); */
/* 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 */