Skip to content
GitLab
Commit b1fff331 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'io_particle_wrapping' into 'master' 

Improvements to i/o and parallel-i/o

See merge request !460
parents 04f8a670 8a35b6f7
Hide whitespace changes
Inline Side-by-side
src/common_io.c
View file @ b1fff331
......@@ -40,10 +40,13 @@
#include "common_io.h"
/* Local includes. */
#include "engine.h"
#include "error.h"
#include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
#include "threadpool.h"
#include "units.h"
#include "version.h"
......@@ -266,7 +269,6 @@ void io_write_attribute_f(hid_t grp, const char* name, float data) {
* @param name The name of the attribute
* @param data The value to write
*/
void io_write_attribute_i(hid_t grp, const char* name, int data) {
io_write_attribute(grp, name, INT, &data, 1);
}
......@@ -295,16 +297,19 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
* @brief Reads the Unit System from an IC file.
* @param h_file The (opened) HDF5 file from which to read.
* @param us The unit_system to fill.
* @param mpi_rank The MPI rank we are on.
*
* If the 'Units' group does not exist in the ICs, cgs units will be assumed
*/
void io_read_unit_system(hid_t h_file, struct unit_system* us) {
void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank) {
/* First check if it exists as this is *not* required. */
const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
if (exists == 0) {
message("'Units' group not found in ICs. Assuming CGS unit system.");
if (mpi_rank == 0)
message("'Units' group not found in ICs. Assuming CGS unit system.");
/* Default to CGS */
us->UnitMass_in_cgs = 1.;
......@@ -319,7 +324,7 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us) {
exists);
}
message("Reading IC units from ICs.");
if (mpi_rank == 0) message("Reading IC units from ICs.");
hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
/* Ok, Read the damn thing */
......@@ -401,6 +406,208 @@ void io_write_code_description(hid_t h_file) {
H5Gclose(h_grpcode);
}
/**
* @brief Mapper function to copy #part or #gpart fields into a buffer.
*/
void io_copy_mapper(void* restrict temp, int N, void* restrict extra_data) {
const struct io_props props = *((const struct io_props*)extra_data);
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
/* How far are we with this chunk? */
char* restrict temp_c = temp;
const ptrdiff_t delta = (temp_c - props.start_temp_c) / copySize;
for (int k = 0; k < N; k++) {
memcpy(&temp_c[k * copySize], props.field + (delta + k) * props.partSize,
copySize);
}
}
/**
* @brief Mapper function to copy #part into a buffer of floats using a
* conversion function.
*/
void io_convert_part_f_mapper(void* restrict temp, int N,
void* restrict extra_data) {
const struct io_props props = *((const struct io_props*)extra_data);
const struct part* restrict parts = props.parts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
/* How far are we with this chunk? */
float* restrict temp_f = temp;
const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
for (int i = 0; i < N; i++)
props.convert_part_f(e, parts + delta + i, &temp_f[i * dim]);
}
/**
* @brief Mapper function to copy #part into a buffer of doubles using a
* conversion function.
*/
void io_convert_part_d_mapper(void* restrict temp, int N,
void* restrict extra_data) {
const struct io_props props = *((const struct io_props*)extra_data);
const struct part* restrict parts = props.parts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
/* How far are we with this chunk? */
double* restrict temp_d = temp;
const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
for (int i = 0; i < N; i++)
props.convert_part_d(e, parts + delta + i, &temp_d[i * dim]);
}
/**
* @brief Mapper function to copy #gpart into a buffer of floats using a
* conversion function.
*/
void io_convert_gpart_f_mapper(void* restrict temp, int N,
void* restrict extra_data) {
const struct io_props props = *((const struct io_props*)extra_data);
const struct gpart* restrict gparts = props.gparts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
/* How far are we with this chunk? */
float* restrict temp_f = temp;
const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
for (int i = 0; i < N; i++)
props.convert_gpart_f(e, gparts + delta + i, &temp_f[i * dim]);
}
/**
* @brief Mapper function to copy #gpart into a buffer of doubles using a
* conversion function.
*/
void io_convert_gpart_d_mapper(void* restrict temp, int N,
void* restrict extra_data) {
const struct io_props props = *((const struct io_props*)extra_data);
const struct gpart* restrict gparts = props.gparts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
/* How far are we with this chunk? */
double* restrict temp_d = temp;
const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
for (int i = 0; i < N; i++)
props.convert_gpart_d(e, gparts + delta + i, &temp_d[i * dim]);
}
/**
* @brief Copy the particle data into a temporary buffer ready for i/o.
*
* @param temp The buffer to be filled. Must be allocated and aligned properly.
* @param e The #engine.
* @param props The #io_props corresponding to the particle field we are
* copying.
* @param N The number of particles to copy
* @param internal_units The system of units used internally.
* @param snapshot_units The system of units used for the snapshots.
*/
void io_copy_temp_buffer(void* temp, const struct engine* e,
struct io_props props, size_t N,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
const size_t num_elements = N * props.dimension;
/* Copy particle data to temporary buffer */
if (props.conversion == 0) { /* No conversion */
/* Prepare some parameters */
char* temp_c = temp;
props.start_temp_c = temp_c;
/* Copy the whole thing into a buffer */
threadpool_map((struct threadpool*)&e->threadpool, io_copy_mapper, temp_c,
N, copySize, 0, (void*)&props);
} else { /* Converting particle to data */
if (props.convert_part_f != NULL) {
/* Prepare some parameters */
float* temp_f = temp;
props.start_temp_f = temp;
props.e = e;
/* Copy the whole thing into a buffer */
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_part_f_mapper, temp_f, N, copySize, 0,
(void*)&props);
} else if (props.convert_part_d != NULL) {
/* Prepare some parameters */
double* temp_d = temp;
props.start_temp_d = temp;
props.e = e;
/* Copy the whole thing into a buffer */
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_part_d_mapper, temp_d, N, copySize, 0,
(void*)&props);
} else if (props.convert_gpart_f != NULL) {
/* Prepare some parameters */
float* temp_f = temp;
props.start_temp_f = temp;
props.e = e;
/* Copy the whole thing into a buffer */
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_gpart_f_mapper, temp_f, N, copySize, 0,
(void*)&props);
} else if (props.convert_gpart_d != NULL) {
/* Prepare some parameters */
double* temp_d = temp;
props.start_temp_d = temp;
props.e = e;
/* Copy the whole thing into a buffer */
threadpool_map((struct threadpool*)&e->threadpool,
io_convert_gpart_d_mapper, temp_d, N, copySize, 0,
(void*)&props);
} else {
error("Missing conversion function");
}
}
/* Unit conversion if necessary */
const double factor =
units_conversion_factor(internal_units, snapshot_units, props.units);
if (factor != 1.) {
/* message("Converting ! factor=%e", factor); */
if (io_is_double_precision(props.type)) {
swift_declare_aligned_ptr(double, temp_d, temp, IO_BUFFER_ALIGNMENT);
for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
} else {
swift_declare_aligned_ptr(float, temp_f, temp, IO_BUFFER_ALIGNMENT);
for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
}
}
}
#endif /* HAVE_HDF5 */
/* ------------------------------------------------------------------------------------------------
......
src/common_io.h
View file @ b1fff331
......@@ -30,6 +30,11 @@
#define FIELD_BUFFER_SIZE 200
#define PARTICLE_GROUP_BUFFER_SIZE 50
#define FILENAME_BUFFER_SIZE 150
#define IO_BUFFER_ALIGNMENT 1024
/* Avoid cyclic inclusion problems */
struct io_props;
struct engine;
#if defined(HAVE_HDF5)
......@@ -68,10 +73,15 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str);
void io_write_code_description(hid_t h_file);
void io_read_unit_system(hid_t h_file, struct unit_system* us);
void io_read_unit_system(hid_t h_file, struct unit_system* us, int mpi_rank);
void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
const char* groupName);
void io_copy_temp_buffer(void* temp, const struct engine* e,
const struct io_props props, size_t N,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units);
#endif /* defined HDF5 */
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
......
src/engine.c
View file @ b1fff331
......@@ -4535,7 +4535,7 @@ void engine_init(struct engine *e, struct space *s,
if (with_aff) engine_unpin();
#endif
if (with_aff) {
if (with_aff && nodeID == 0) {
#ifdef HAVE_SETAFFINITY
#ifdef WITH_MPI
printf("[%04i] %s engine_init: cpu map is [ ", nodeID,
......
src/gravity/Default/gravity_io.h
View file @ b1fff331
......@@ -21,6 +21,20 @@
#include "io_properties.h"
void convert_gpart_pos(const struct engine* e, const struct gpart* gp,
double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(gp->x[0], 0.0, e->s->dim[0]);
ret[1] = box_wrap(gp->x[1], 0.0, e->s->dim[1]);
ret[2] = box_wrap(gp->x[2], 0.0, e->s->dim[2]);
} else {
ret[0] = gp->x[0];
ret[1] = gp->x[1];
ret[2] = gp->x[2];
}
}
/**
* @brief Specifies which g-particle fields to read from a dataset
*
......@@ -52,15 +66,15 @@ void darkmatter_read_particles(struct gpart* gparts, struct io_props* list,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void darkmatter_write_particles(struct gpart* gparts, struct io_props* list,
int* num_fields) {
void darkmatter_write_particles(const struct gpart* gparts,
struct io_props* list, int* num_fields) {
/* Say how much we want to read */
*num_fields = 5;
/* List what we want to read */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
gparts, x);
list[0] = io_make_output_field_convert_gpart(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
gparts, v_full);
list[2] =
......
src/hydro/Default/hydro_io.h
View file @ b1fff331
......@@ -53,6 +53,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
} else {
ret[0] = p->x[0];
ret[1] = p->x[1];
ret[2] = p->x[2];
}
}
/**
* @brief Specifies which particle fields to write to a dataset
*
......@@ -66,8 +80,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 8;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
parts, x);
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[2] =
......
src/hydro/Default/hydro_part.h
View file @ b1fff331
......@@ -125,6 +125,16 @@ struct part {
/* Particle time-bin */
timebin_t time_bin;
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
integertime_t ti_drift;
/* Time of the last kick */
integertime_t ti_kick;
#endif
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_HYDRO_PART_H */
src/hydro/Gadget2/hydro_io.h
View file @ b1fff331
......@@ -55,14 +55,28 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
float convert_u(struct engine* e, struct part* p) {
void convert_u(const struct engine* e, const struct part* p, float* ret) {
return hydro_get_internal_energy(p);
ret[0] = hydro_get_internal_energy(p);
}
float convert_P(struct engine* e, struct part* p) {
void convert_P(const struct engine* e, const struct part* p, float* ret) {
return hydro_get_pressure(p);
ret[0] = hydro_get_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
} else {
ret[0] = p->x[0];
ret[1] = p->x[1];
ret[2] = p->x[2];
}
}
/**
......@@ -72,14 +86,14 @@ float convert_P(struct engine* e, struct part* p) {
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
void hydro_write_particles(struct part* parts, struct io_props* list,
void hydro_write_particles(const struct part* parts, struct io_props* list,
int* num_fields) {
*num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
parts, x);
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[2] =
......@@ -96,9 +110,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[8] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, rho, convert_u);
parts, convert_u);
list[9] = io_make_output_field_convert_part(
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
"Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
}
/**
......
src/hydro/Gizmo/hydro_io.h
View file @ b1fff331
......@@ -69,10 +69,11 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
*
* @param e #engine.
* @param p Particle.
* @return Internal energy of the particle
* @param ret (return) Internal energy of the particle
*/
float convert_u(struct engine* e, struct part* p) {
return hydro_get_internal_energy(p);
void convert_u(const struct engine* e, const struct part* p, float* ret) {
ret[0] = hydro_get_internal_energy(p);
}
/**
......@@ -80,10 +81,10 @@ float convert_u(struct engine* e, struct part* p) {
*
* @param e #engine.
* @param p Particle.
* @return Entropic function of the particle
* @param ret (return) Entropic function of the particle
*/
float convert_A(struct engine* e, struct part* p) {
return hydro_get_entropy(p);
void convert_A(const struct engine* e, const struct part* p, float* ret) {
ret[0] = hydro_get_entropy(p);
}
/**
......@@ -93,9 +94,9 @@ float convert_A(struct engine* e, struct part* p) {
* @param p Particle.
* @return Total energy of the particle
*/
float convert_Etot(struct engine* e, struct part* p) {
void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
#ifdef GIZMO_TOTAL_ENERGY
return p->conserved.energy;
ret[0] = p->conserved.energy;
#else
float momentum2;
......@@ -103,10 +104,24 @@ float convert_Etot(struct engine* e, struct part* p) {
p->conserved.momentum[1] * p->conserved.momentum[1] +
p->conserved.momentum[2] * p->conserved.momentum[2];
return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
ret[0] = p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
#endif
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
} else {
ret[0] = p->x[0];
ret[1] = p->x[1];
ret[2] = p->x[2];
}
}
/**
* @brief Specifies which particle fields to write to a dataset
*
......@@ -120,8 +135,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
parts, x);
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
primitives.v);
list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
......@@ -130,18 +145,17 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
parts, h);
list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
parts, primitives.P, convert_u);
parts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
primitives.rho);
list[7] = io_make_output_field_convert_part(
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, primitives.P, convert_A);
"Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, convert_A);
list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, primitives.P);
list[9] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
list[9] = io_make_output_field_convert_part(
"TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, convert_Etot);
}
/**
......
src/hydro/Minimal/hydro_io.h
View file @ b1fff331
......@@ -69,14 +69,28 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
float convert_S(struct engine* e, struct part* p) {
void convert_S(const struct engine* e, const struct part* p, float* ret) {
return hydro_get_entropy(p);
ret[0] = hydro_get_entropy(p);
}
float convert_P(struct engine* e, struct part* p) {
void convert_P(const struct engine* e, const struct part* p, float* ret) {
return hydro_get_pressure(p);
ret[0] = hydro_get_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
ret[1] = box_wrap(p->x[1], 0.0, e->s->dim[1]);
ret[2] = box_wrap(p->x[2], 0.0, e->s->dim[2]);
} else {
ret[0] = p->x[0];
ret[1] = p->x[1];
ret[2] = p->x[2];
}
}
/**
......@@ -92,8 +106,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 10;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
parts, x);
list[0] = io_make_output_field_convert_part(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[2] =
......@@ -108,11 +122,10 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_ACCELERATION, parts, a_hydro);
list[7] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[8] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
UNIT_CONV_ENTROPY_PER_UNIT_MASS,