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

Merge branch 'master' into gravity_fixes

parents 55f708b6 cb0b4486
...@@ -349,10 +349,12 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted( ...@@ -349,10 +349,12 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
y[i] = (float)(parts_i[idx].x[1] - total_ci_shift[1]); y[i] = (float)(parts_i[idx].x[1] - total_ci_shift[1]);
z[i] = (float)(parts_i[idx].x[2] - total_ci_shift[2]); z[i] = (float)(parts_i[idx].x[2] - total_ci_shift[2]);
h[i] = parts_i[idx].h; h[i] = parts_i[idx].h;
m[i] = parts_i[idx].mass;
vx[i] = parts_i[idx].v[0]; vx[i] = parts_i[idx].v[0];
vy[i] = parts_i[idx].v[1]; vy[i] = parts_i[idx].v[1];
vz[i] = parts_i[idx].v[2]; vz[i] = parts_i[idx].v[2];
#ifdef GADGET2_SPH
m[i] = parts_i[idx].mass;
#endif
} }
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
...@@ -431,10 +433,12 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted( ...@@ -431,10 +433,12 @@ __attribute__((always_inline)) INLINE void cache_read_two_partial_cells_sorted(
yj[i] = (float)(parts_j[idx].x[1] - total_cj_shift[1]); yj[i] = (float)(parts_j[idx].x[1] - total_cj_shift[1]);
zj[i] = (float)(parts_j[idx].x[2] - total_cj_shift[2]); zj[i] = (float)(parts_j[idx].x[2] - total_cj_shift[2]);
hj[i] = parts_j[idx].h; hj[i] = parts_j[idx].h;
mj[i] = parts_j[idx].mass;
vxj[i] = parts_j[idx].v[0]; vxj[i] = parts_j[idx].v[0];
vyj[i] = parts_j[idx].v[1]; vyj[i] = parts_j[idx].v[1];
vzj[i] = parts_j[idx].v[2]; vzj[i] = parts_j[idx].v[2];
#ifdef GADGET2_SPH
mj[i] = parts_j[idx].mass;
#endif
} }
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
...@@ -572,15 +576,17 @@ cache_read_two_partial_cells_sorted_force( ...@@ -572,15 +576,17 @@ cache_read_two_partial_cells_sorted_force(
y[i] = (float)(parts_i[idx].x[1] - total_ci_shift[1]); y[i] = (float)(parts_i[idx].x[1] - total_ci_shift[1]);
z[i] = (float)(parts_i[idx].x[2] - total_ci_shift[2]); z[i] = (float)(parts_i[idx].x[2] - total_ci_shift[2]);
h[i] = parts_i[idx].h; h[i] = parts_i[idx].h;
m[i] = parts_i[idx].mass;
vx[i] = parts_i[idx].v[0]; vx[i] = parts_i[idx].v[0];
vy[i] = parts_i[idx].v[1]; vy[i] = parts_i[idx].v[1];
vz[i] = parts_i[idx].v[2]; vz[i] = parts_i[idx].v[2];
#ifdef GADGET2_SPH
m[i] = parts_i[idx].mass;
rho[i] = parts_i[idx].rho; rho[i] = parts_i[idx].rho;
grad_h[i] = parts_i[idx].force.f; grad_h[i] = parts_i[idx].force.f;
pOrho2[i] = parts_i[idx].force.P_over_rho2; pOrho2[i] = parts_i[idx].force.P_over_rho2;
balsara[i] = parts_i[idx].force.balsara; balsara[i] = parts_i[idx].force.balsara;
soundspeed[i] = parts_i[idx].force.soundspeed; soundspeed[i] = parts_i[idx].force.soundspeed;
#endif
} }
/* Pad cache with fake particles that exist outside the cell so will not /* Pad cache with fake particles that exist outside the cell so will not
...@@ -635,15 +641,17 @@ cache_read_two_partial_cells_sorted_force( ...@@ -635,15 +641,17 @@ cache_read_two_partial_cells_sorted_force(
yj[i] = (float)(parts_j[idx].x[1] - total_cj_shift[1]); yj[i] = (float)(parts_j[idx].x[1] - total_cj_shift[1]);
zj[i] = (float)(parts_j[idx].x[2] - total_cj_shift[2]); zj[i] = (float)(parts_j[idx].x[2] - total_cj_shift[2]);
hj[i] = parts_j[idx].h; hj[i] = parts_j[idx].h;
mj[i] = parts_j[idx].mass;
vxj[i] = parts_j[idx].v[0]; vxj[i] = parts_j[idx].v[0];
vyj[i] = parts_j[idx].v[1]; vyj[i] = parts_j[idx].v[1];
vzj[i] = parts_j[idx].v[2]; vzj[i] = parts_j[idx].v[2];
#ifdef GADGET2_SPH
mj[i] = parts_j[idx].mass;
rhoj[i] = parts_j[idx].rho; rhoj[i] = parts_j[idx].rho;
grad_hj[i] = parts_j[idx].force.f; grad_hj[i] = parts_j[idx].force.f;
pOrho2j[i] = parts_j[idx].force.P_over_rho2; pOrho2j[i] = parts_j[idx].force.P_over_rho2;
balsaraj[i] = parts_j[idx].force.balsara; balsaraj[i] = parts_j[idx].force.balsara;
soundspeedj[i] = parts_j[idx].force.soundspeed; soundspeedj[i] = parts_j[idx].force.soundspeed;
#endif
} }
/* Pad cache with fake particles that exist outside the cell so will not /* Pad cache with fake particles that exist outside the cell so will not
......
...@@ -40,10 +40,13 @@ ...@@ -40,10 +40,13 @@
#include "common_io.h" #include "common_io.h"
/* Local includes. */ /* Local includes. */
#include "engine.h"
#include "error.h" #include "error.h"
#include "hydro.h" #include "hydro.h"
#include "io_properties.h"
#include "kernel_hydro.h" #include "kernel_hydro.h"
#include "part.h" #include "part.h"
#include "threadpool.h"
#include "units.h" #include "units.h"
#include "version.h" #include "version.h"
...@@ -266,7 +269,6 @@ void io_write_attribute_f(hid_t grp, const char* name, float data) { ...@@ -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 name The name of the attribute
* @param data The value to write * @param data The value to write
*/ */
void io_write_attribute_i(hid_t grp, const char* name, int data) { void io_write_attribute_i(hid_t grp, const char* name, int data) {
io_write_attribute(grp, name, INT, &data, 1); io_write_attribute(grp, name, INT, &data, 1);
} }
...@@ -295,15 +297,18 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) { ...@@ -295,15 +297,18 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str) {
* @brief Reads the Unit System from an IC file. * @brief Reads the Unit System from an IC file.
* @param h_file The (opened) HDF5 file from which to read. * @param h_file The (opened) HDF5 file from which to read.
* @param us The unit_system to fill. * @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 * 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. */ /* First check if it exists as this is *not* required. */
const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT); const htri_t exists = H5Lexists(h_file, "/Units", H5P_DEFAULT);
if (exists == 0) { if (exists == 0) {
if (mpi_rank == 0)
message("'Units' group not found in ICs. Assuming CGS unit system."); message("'Units' group not found in ICs. Assuming CGS unit system.");
/* Default to CGS */ /* Default to CGS */
...@@ -319,7 +324,7 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us) { ...@@ -319,7 +324,7 @@ void io_read_unit_system(hid_t h_file, struct unit_system* us) {
exists); 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); hid_t h_grp = H5Gopen(h_file, "/Units", H5P_DEFAULT);
/* Ok, Read the damn thing */ /* Ok, Read the damn thing */
...@@ -401,6 +406,208 @@ void io_write_code_description(hid_t h_file) { ...@@ -401,6 +406,208 @@ void io_write_code_description(hid_t h_file) {
H5Gclose(h_grpcode); 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 */ #endif /* HAVE_HDF5 */
/* ------------------------------------------------------------------------------------------------ /* ------------------------------------------------------------------------------------------------
......
...@@ -30,6 +30,11 @@ ...@@ -30,6 +30,11 @@
#define FIELD_BUFFER_SIZE 200 #define FIELD_BUFFER_SIZE 200
#define PARTICLE_GROUP_BUFFER_SIZE 50 #define PARTICLE_GROUP_BUFFER_SIZE 50
#define FILENAME_BUFFER_SIZE 150 #define FILENAME_BUFFER_SIZE 150
#define IO_BUFFER_ALIGNMENT 1024
/* Avoid cyclic inclusion problems */
struct io_props;
struct engine;
#if defined(HAVE_HDF5) #if defined(HAVE_HDF5)
...@@ -68,10 +73,15 @@ void io_write_attribute_s(hid_t grp, const char* name, const char* str); ...@@ -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_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, void io_write_unit_system(hid_t h_grp, const struct unit_system* us,
const char* groupName); 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 */ #endif /* defined HDF5 */
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot, void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
......
...@@ -2624,9 +2624,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements, ...@@ -2624,9 +2624,8 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Self-interaction? */ /* Self-interaction? */
else if (t->type == task_type_self && t->subtype == task_subtype_density) { else if (t->type == task_type_self && t->subtype == task_subtype_density) {
/* Make all density tasks depend on the drift and the sorts. */ /* Make the self-density tasks depend on the drift only. */
scheduler_addunlock(sched, t->ci->super->drift_part, t); scheduler_addunlock(sched, t->ci->super->drift_part, t);
scheduler_addunlock(sched, t->ci->super->sorts, t);
#ifdef EXTRA_HYDRO_LOOP #ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop. */ /* Start by constructing the task for the second and third hydro loop. */
...@@ -5004,7 +5003,7 @@ void engine_init(struct engine *e, struct space *s, ...@@ -5004,7 +5003,7 @@ void engine_init(struct engine *e, struct space *s,
if (with_aff) engine_unpin(); if (with_aff) engine_unpin();
#endif #endif
if (with_aff) { if (with_aff && nodeID == 0) {
#ifdef HAVE_SETAFFINITY #ifdef HAVE_SETAFFINITY
#ifdef WITH_MPI #ifdef WITH_MPI
printf("[%04i] %s engine_init: cpu map is [ ", nodeID, printf("[%04i] %s engine_init: cpu map is [ ", nodeID,
......
...@@ -21,6 +21,20 @@ ...@@ -21,6 +21,20 @@
#include "io_properties.h" #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 * @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, ...@@ -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 list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write. * @param num_fields The number of i/o fields to write.
*/ */
void darkmatter_write_particles(struct gpart* gparts, struct io_props* list, void darkmatter_write_particles(const struct gpart* gparts,
int* num_fields) { struct io_props* list, int* num_fields) {
/* Say how much we want to read */ /* Say how much we want to read */
*num_fields = 5; *num_fields = 5;
/* List what we want to read */ /* List what we want to read */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, list[0] = io_make_output_field_convert_gpart(
gparts, x); "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
gparts, v_full); gparts, v_full);
list[2] = list[2] =
......
...@@ -53,6 +53,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list, ...@@ -53,6 +53,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho); 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 * @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, ...@@ -66,8 +80,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
*num_fields = 8; *num_fields = 8;
/* List what we want to write */ /* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, list[0] = io_make_output_field_convert_part(
parts, x); "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
list[1] = list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v); io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
list[2] = list[2] =
......
...@@ -125,6 +125,16 @@ struct part { ...@@ -125,6 +125,16 @@ struct part {
/* Particle time-bin */ /* Particle time-bin */
timebin_t 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; } SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_HYDRO_PART_H */ #endif /* SWIFT_DEFAULT_HYDRO_PART_H */
...@@ -55,14 +55,28 @@ void hydro_read_particles(struct part* parts, struct io_props* list, ...@@ -55,14 +55,28 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho); 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,