diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h index 26eed8a0a680e13130dd257d8b31e6cd00392f6d..671ccc24e95f322d71a734c0a39eface2a93bcec 100644 --- a/src/gravity/Default/gravity_io.h +++ b/src/gravity/Default/gravity_io.h @@ -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] = diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index 3e46b2351eb6a3871946dd9c69c8d108b10da0df..d89fe5c277b0273e53304d77a241e14d2c71fc5d 100644 --- a/src/hydro/Gadget2/hydro_io.h +++ b/src/hydro/Gadget2/hydro_io.h @@ -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); } /** diff --git a/src/io_properties.h b/src/io_properties.h index 9fcf1a1ac67cae6afab6870369e51d06c752fc11..ff8ad7a32e7e4cbb75f6a1f4ea2ad3d714a64be6 100644 --- a/src/io_properties.h +++ b/src/io_properties.h @@ -16,7 +16,6 @@ * along with this program. If not, see <http://www.gnu.org/licenses/>. * ******************************************************************************/ - #ifndef SWIFT_IO_PROPERTIES_H #define SWIFT_IO_PROPERTIES_H @@ -29,6 +28,16 @@ */ enum DATA_IMPORTANCE { COMPULSORY = 1, OPTIONAL = 0, UNUSED = -1 }; +/* Helper typedefs */ +typedef void (*conversion_func_part_float)(const struct engine*, + const struct part*, float*); +typedef void (*conversion_func_part_double)(const struct engine*, + const struct part*, double*); +typedef void (*conversion_func_gpart_float)(const struct engine*, + const struct gpart*, float*); +typedef void (*conversion_func_gpart_double)(const struct engine*, + const struct gpart*, double*); + /** * @brief The properties of a given dataset for i/o */ @@ -56,14 +65,19 @@ struct io_props { size_t partSize; /* The particle arrays */ - struct part* parts; - struct gpart* gparts; + const struct part* parts; + const struct gpart* gparts; + + /* Are we converting? */ + int conversion; /* Conversion function for part */ - float (*convert_part)(struct engine*, struct part*); + conversion_func_part_float convert_part_f; + conversion_func_part_double convert_part_d; /* Conversion function for gpart */ - float (*convert_gpart)(struct engine*, struct gpart*); + conversion_func_gpart_float convert_gpart_f; + conversion_func_gpart_double convert_gpart_d; }; /** @@ -101,8 +115,11 @@ struct io_props io_make_input_field_(char name[FIELD_BUFFER_SIZE], r.partSize = partSize; r.parts = NULL; r.gparts = NULL; - r.convert_part = NULL; - r.convert_gpart = NULL; + r.conversion = 0; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; return r; } @@ -140,8 +157,11 @@ struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE], r.partSize = partSize; r.parts = NULL; r.gparts = NULL; - r.convert_part = NULL; - r.convert_gpart = NULL; + r.conversion = 0; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; return r; } @@ -149,11 +169,10 @@ struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE], /** * @brief Constructs an #io_props (with conversion) from its parameters */ -#define io_make_output_field_convert_part(name, type, dim, units, part, field, \ - convert) \ - io_make_output_field_convert_part_(name, type, dim, units, \ - (char*)(&(part[0]).field), \ - sizeof(part[0]), part, convert) +#define io_make_output_field_convert_part(name, type, dim, units, part, \ + convert) \ + io_make_output_field_convert_part_##type(name, type, dim, units, \ + sizeof(part[0]), part, convert) /** * @brief Construct an #io_props from its parameters @@ -162,17 +181,16 @@ struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE], * @param type The type of the data * @param dimension Dataset dimension (1D, 3D, ...) * @param units The units of the dataset - * @param field Pointer to the field of the first particle * @param partSize The size in byte of the particle * @param parts The particle array * @param functionPtr The function used to convert a particle to a float * * Do not call this function directly. Use the macro defined above. */ -struct io_props io_make_output_field_convert_part_( +struct io_props io_make_output_field_convert_part_FLOAT( char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, - enum unit_conversion_factor units, char* field, size_t partSize, - struct part* parts, float (*functionPtr)(struct engine*, struct part*)) { + enum unit_conversion_factor units, size_t partSize, + const struct part* parts, conversion_func_part_float functionPtr) { struct io_props r; strcpy(r.name, name); @@ -180,12 +198,52 @@ struct io_props io_make_output_field_convert_part_( r.dimension = dimension; r.importance = UNUSED; r.units = units; - r.field = field; + r.field = NULL; r.partSize = partSize; r.parts = parts; r.gparts = NULL; - r.convert_part = functionPtr; - r.convert_gpart = NULL; + r.conversion = 1; + r.convert_part_f = functionPtr; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; + + return r; +} + +/** + * @brief Construct an #io_props from its parameters + * + * @param name Name of the field to read + * @param type The type of the data + * @param dimension Dataset dimension (1D, 3D, ...) + * @param units The units of the dataset + * @param partSize The size in byte of the particle + * @param parts The particle array + * @param functionPtr The function used to convert a particle to a float + * + * Do not call this function directly. Use the macro defined above. + */ +struct io_props io_make_output_field_convert_part_DOUBLE( + char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, + enum unit_conversion_factor units, size_t partSize, + const struct part* parts, conversion_func_part_double functionPtr) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = UNUSED; + r.units = units; + r.field = NULL; + r.partSize = partSize; + r.parts = parts; + r.gparts = NULL; + r.conversion = 1; + r.convert_part_f = NULL; + r.convert_part_d = functionPtr; + r.convert_gpart_f = NULL; + r.convert_gpart_d = NULL; return r; } @@ -193,11 +251,10 @@ struct io_props io_make_output_field_convert_part_( /** * @brief Constructs an #io_props (with conversion) from its parameters */ -#define io_make_output_field_convert_gpart(name, type, dim, units, part, \ - field, convert) \ - io_make_output_field_convert_gpart_(name, type, dim, units, \ - (char*)(&(part[0]).field), \ - sizeof(part[0]), gpart, convert) +#define io_make_output_field_convert_gpart(name, type, dim, units, gpart, \ + convert) \ + io_make_output_field_convert_gpart_##type(name, type, dim, units, \ + sizeof(gpart[0]), gpart, convert) /** * @brief Construct an #io_props from its parameters @@ -206,17 +263,16 @@ struct io_props io_make_output_field_convert_part_( * @param type The type of the data * @param dimension Dataset dimension (1D, 3D, ...) * @param units The units of the dataset - * @param field Pointer to the field of the first particle - * @param partSize The size in byte of the particle + * @param gpartSize The size in byte of the particle * @param gparts The particle array * @param functionPtr The function used to convert a g-particle to a float * * Do not call this function directly. Use the macro defined above. */ -struct io_props io_make_output_field_convert_gpart_( +struct io_props io_make_output_field_convert_gpart_FLOAT( char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, - enum unit_conversion_factor units, char* field, size_t partSize, - struct gpart* gparts, float (*functionPtr)(struct engine*, struct gpart*)) { + enum unit_conversion_factor units, size_t gpartSize, + const struct gpart* gparts, conversion_func_gpart_float functionPtr) { struct io_props r; strcpy(r.name, name); @@ -224,12 +280,52 @@ struct io_props io_make_output_field_convert_gpart_( r.dimension = dimension; r.importance = UNUSED; r.units = units; - r.field = field; - r.partSize = partSize; + r.field = NULL; + r.partSize = gpartSize; + r.parts = NULL; + r.gparts = gparts; + r.conversion = 1; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = functionPtr; + r.convert_gpart_d = NULL; + + return r; +} + +/** + * @brief Construct an #io_props from its parameters + * + * @param name Name of the field to read + * @param type The type of the data + * @param dimension Dataset dimension (1D, 3D, ...) + * @param units The units of the dataset + * @param gpartSize The size in byte of the particle + * @param gparts The particle array + * @param functionPtr The function used to convert a g-particle to a float + * + * Do not call this function directly. Use the macro defined above. + */ +struct io_props io_make_output_field_convert_gpart_DOUBLE( + char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension, + enum unit_conversion_factor units, size_t gpartSize, + const struct gpart* gparts, conversion_func_gpart_double functionPtr) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = UNUSED; + r.units = units; + r.field = NULL; + r.partSize = gpartSize; r.parts = NULL; r.gparts = gparts; - r.convert_part = NULL; - r.convert_gpart = functionPtr; + r.conversion = 1; + r.convert_part_f = NULL; + r.convert_part_d = NULL; + r.convert_gpart_f = NULL; + r.convert_gpart_d = functionPtr; return r; } diff --git a/src/serial_io.c b/src/serial_io.c index eb1e0e23fb34fd8d6a21230d9e38cfe82c47df1d..b161a2f44b34308b3e0cb5c4fcaa55d884b896d6 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -176,9 +176,9 @@ void readArray(hid_t grp, const struct io_props props, size_t N, * Routines writing an output file *-----------------------------------------------------------------------------*/ -void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, - char* partTypeGroupName, const struct io_props props, - unsigned long long N_total, +void prepareArray(const struct engine* e, hid_t grp, char* fileName, + FILE* xmfFile, char* partTypeGroupName, + const struct io_props props, unsigned long long N_total, const struct unit_system* internal_units, const struct unit_system* snapshot_units) { @@ -282,15 +282,17 @@ void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, * @todo A better version using HDF5 hyper-slabs to write the file directly from * the part array will be written once the structures have been stabilized. */ -void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, - char* partTypeGroupName, const struct io_props props, size_t N, - long long N_total, int mpi_rank, long long offset, +void writeArray(const struct engine* e, hid_t grp, char* fileName, + FILE* xmfFile, char* partTypeGroupName, + const struct io_props props, size_t N, long long N_total, + int mpi_rank, long long offset, 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; + const size_t dim = props.dimension; /* message("Writing '%s' array...", props.name); */ @@ -300,28 +302,63 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, internal_units, snapshot_units); /* Allocate temporary buffer */ - void* temp = malloc(num_elements * io_sizeof_type(props.type)); - if (temp == NULL) error("Unable to allocate memory for temporary buffer"); + void* temp = NULL; + if (posix_memalign((void**)&temp, SWIFT_CACHE_ALIGNMENT, + num_elements * io_sizeof_type(props.type)) != 0) + error("Unable to allocate temporary i/o buffer"); /* Copy particle data to temporary buffer */ - if (props.convert_part == NULL && - props.convert_gpart == NULL) { /* No conversion */ + if (props.conversion == 0) { /* No conversion */ char* temp_c = temp; for (size_t i = 0; i < N; ++i) memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize); - } else if (props.convert_part != NULL) { /* conversion (for parts)*/ + } else { /* Converting particle to data */ - float* temp_f = temp; - for (size_t i = 0; i < N; ++i) - temp_f[i] = props.convert_part(e, &props.parts[i]); + if (props.convert_part_f != NULL) { - } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/ + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct part, parts, props.parts, + SWIFT_STRUCT_ALIGNMENT); - float* temp_f = temp; - for (size_t i = 0; i < N; ++i) - temp_f[i] = props.convert_gpart(e, &props.gparts[i]); + /* float conversion for parts */ + for (size_t i = 0; i < N; i++) + props.convert_part_f(e, &parts[i], &temp_f[i * dim]); + + } else if (props.convert_part_d != NULL) { + + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct part, parts, props.parts, + SWIFT_STRUCT_ALIGNMENT); + + /* double conversion for parts */ + for (size_t i = 0; i < N; i++) + props.convert_part_d(e, &parts[i], &temp_d[i * dim]); + + } else if (props.convert_gpart_f != NULL) { + + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct gpart, gparts, props.gparts, + SWIFT_STRUCT_ALIGNMENT); + + /* float conversion for gparts */ + for (size_t i = 0; i < N; i++) + props.convert_gpart_f(e, &gparts[i], &temp_f[i * dim]); + + } else if (props.convert_gpart_d != NULL) { + + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct gpart, gparts, props.gparts, + SWIFT_STRUCT_ALIGNMENT); + + /* double conversion for gparts */ + for (size_t i = 0; i < N; i++) + props.convert_gpart_d(e, &gparts[i], &temp_d[i * dim]); + + } else { + error("Missing conversion function"); + } } /* Unit conversion if necessary */ @@ -332,10 +369,10 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, /* message("Converting ! factor=%e", factor); */ if (io_is_double_precision(props.type)) { - double* temp_d = temp; + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor; } else { - float* temp_f = temp; + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor; } } diff --git a/src/single_io.c b/src/single_io.c index 194563352dff5570b8703f828fac95bccbf7409f..e78293ba317149d469ef8d47e04279d9f076256b 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -169,40 +169,77 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N, * @todo A better version using HDF5 hyper-slabs to write the file directly from * the part array will be written once the structures have been stabilized. */ -void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, - char* partTypeGroupName, const struct io_props props, size_t N, +void writeArray(const struct engine* e, hid_t grp, char* fileName, + FILE* xmfFile, char* partTypeGroupName, + const 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; + const size_t dim = props.dimension; /* message("Writing '%s' array...", props.name); */ /* Allocate temporary buffer */ - void* temp = malloc(num_elements * io_sizeof_type(props.type)); - if (temp == NULL) error("Unable to allocate memory for temporary buffer"); + void* temp = NULL; + if (posix_memalign((void**)&temp, SWIFT_CACHE_ALIGNMENT, + num_elements * io_sizeof_type(props.type)) != 0) + error("Unable to allocate temporary i/o buffer"); /* Copy particle data to temporary buffer */ - if (props.convert_part == NULL && - props.convert_gpart == NULL) { /* No conversion */ + if (props.conversion == 0) { /* No conversion */ char* temp_c = temp; for (size_t i = 0; i < N; ++i) memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize); - } else if (props.convert_part != NULL) { /* conversion (for parts)*/ + } else { /* Converting particle to data */ - float* temp_f = temp; - for (size_t i = 0; i < N; ++i) - temp_f[i] = props.convert_part(e, &props.parts[i]); + if (props.convert_part_f != NULL) { - } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/ + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct part, parts, props.parts, + SWIFT_STRUCT_ALIGNMENT); - float* temp_f = temp; - for (size_t i = 0; i < N; ++i) - temp_f[i] = props.convert_gpart(e, &props.gparts[i]); + /* float conversion for parts */ + for (size_t i = 0; i < N; i++) + props.convert_part_f(e, &parts[i], &temp_f[i * dim]); + + } else if (props.convert_part_d != NULL) { + + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct part, parts, props.parts, + SWIFT_STRUCT_ALIGNMENT); + + /* double conversion for parts */ + for (size_t i = 0; i < N; i++) + props.convert_part_d(e, &parts[i], &temp_d[i * dim]); + + } else if (props.convert_gpart_f != NULL) { + + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct gpart, gparts, props.gparts, + SWIFT_STRUCT_ALIGNMENT); + + /* float conversion for gparts */ + for (size_t i = 0; i < N; i++) + props.convert_gpart_f(e, &gparts[i], &temp_f[i * dim]); + + } else if (props.convert_gpart_d != NULL) { + + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); + swift_declare_aligned_ptr(const struct gpart, gparts, props.gparts, + SWIFT_STRUCT_ALIGNMENT); + + /* double conversion for gparts */ + for (size_t i = 0; i < N; i++) + props.convert_gpart_d(e, &gparts[i], &temp_d[i * dim]); + + } else { + error("Missing conversion function"); + } } /* Unit conversion if necessary */ @@ -213,10 +250,10 @@ void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, /* message("Converting ! factor=%e", factor); */ if (io_is_double_precision(props.type)) { - double* temp_d = temp; + swift_declare_aligned_ptr(double, temp_d, temp, SWIFT_CACHE_ALIGNMENT); for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor; } else { - float* temp_f = temp; + swift_declare_aligned_ptr(float, temp_f, temp, SWIFT_CACHE_ALIGNMENT); for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor; } }