diff --git a/src/cell.c b/src/cell.c index cd83ddab570a3c848e0502629534743e8f383a4f..3fcb8c45633926d71998f1ebbbab6f4162667b2e 100644 --- a/src/cell.c +++ b/src/cell.c @@ -638,11 +638,11 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) { */ void cell_init_parts(struct cell *c, void *data) { - struct part *p = c->parts; - struct xpart *xp = c->xparts; - const int count = c->count; + struct part *restrict p = c->parts; + struct xpart *restrict xp = c->xparts; + const size_t count = c->count; - for (int i = 0; i < count; ++i) { + for (size_t i = 0; i < count; ++i) { p[i].ti_begin = 0; p[i].ti_end = 0; xp[i].v_full[0] = p[i].v[0]; @@ -665,13 +665,14 @@ void cell_init_parts(struct cell *c, void *data) { */ void cell_init_gparts(struct cell *c, void *data) { - struct gpart *gp = c->gparts; - const int gcount = c->gcount; + struct gpart *restrict gp = c->gparts; + const size_t gcount = c->gcount; - for (int i = 0; i < gcount; ++i) { + for (size_t i = 0; i < gcount; ++i) { gp[i].ti_begin = 0; gp[i].ti_end = 0; gravity_first_init_gpart(&gp[i]); + gravity_init_gpart(&gp[i]); } c->ti_end_min = 0; c->ti_end_max = 0; diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h index 0f62511eced181fbf3b5b781a50314dd08f7c0ef..7ef6db76e47f5a1cb2efee9d90fabb9f2ca7dd7f 100644 --- a/src/gravity/Default/gravity.h +++ b/src/gravity/Default/gravity.h @@ -82,7 +82,7 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart( * * @param gp The particle to act upon */ -__attribute__((always_inline)) INLINE static void gravity_init_part( +__attribute__((always_inline)) INLINE static void gravity_init_gpart( struct gpart* gp) { /* Zero the acceleration */ diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index ee693d60999779d5a3cd888e4737208868879ba8..a316eb21db364b255ce9ccb4f97f39eb61e0330f 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -136,15 +136,12 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Some smoothing length multiples. */ const float h = p->h; const float ih = 1.0f / h; - const float ih2 = ih * ih; - const float ih4 = ih2 * ih2; /* Pre-compute some stuff for the balsara switch. */ - const float normDiv_v = fabs(p->density.div_v / p->rho * ih4); + const float normDiv_v = fabs(p->density.div_v); const float normRot_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + p->density.rot_v[1] * p->density.rot_v[1] + - p->density.rot_v[2] * p->density.rot_v[2]) / - p->rho * ih4; + p->density.rot_v[2] * p->density.rot_v[2]); /* Compute this particle's sound speed. */ const float u = p->u; diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h index a122fc24961644073470f9638bad967a52ef4926..1c3e70753d417eb59d0fda5a7acfd1c24ab93045 100644 --- a/src/hydro/Default/hydro_io.h +++ b/src/hydro/Default/hydro_io.h @@ -42,7 +42,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list, list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, UNIT_CONV_LENGTH, parts, h); list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, - UNIT_CONV_ENERGY, parts, u); + UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, parts, id); list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL, @@ -72,7 +72,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list, io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass); list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, parts, h); - list[4] = io_make_output_field("Entropy", FLOAT, 1, + list[4] = io_make_output_field("InternalEnergy", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, parts, id); diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h index 946145587262a1502f9e2ee847f10500c1a6b986..fa207a7fb9ded16d6618df0cb448d1dd49088372 100644 --- a/src/hydro/Default/hydro_part.h +++ b/src/hydro/Default/hydro_part.h @@ -69,43 +69,43 @@ struct part { float alpha; /* Store density/force specific stuff. */ - // union { + union { - struct { + struct { - /* Particle velocity divergence. */ - float div_v; + /* Particle velocity divergence. */ + float div_v; - /* Derivative of particle number density. */ - float wcount_dh; + /* Derivative of particle number density. */ + float wcount_dh; - /* Particle velocity curl. */ - float rot_v[3]; + /* Particle velocity curl. */ + float rot_v[3]; - /* Particle number density. */ - float wcount; + /* Particle number density. */ + float wcount; - } density; + } density; - struct { + struct { - /* Balsara switch */ - float balsara; + /* Balsara switch */ + float balsara; - /* Aggregate quantities. */ - float POrho2; + /* Aggregate quantities. */ + float POrho2; - /* Change in particle energy over time. */ - float u_dt; + /* Change in particle energy over time. */ + float u_dt; - /* Signal velocity */ - float v_sig; + /* Signal velocity */ + float v_sig; - /* Sound speed */ - float c; + /* Sound speed */ + float c; - } force; - //}; + } force; + }; /* Particle mass. */ float mass; diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h index c8db0b3b64cb03e4e013c6dd8daeb16a138b6709..1614a483ac3007d37770e04d40833f78d19a8cd2 100644 --- a/src/hydro/Gadget2/hydro_io.h +++ b/src/hydro/Gadget2/hydro_io.h @@ -17,6 +17,7 @@ * ******************************************************************************/ +#include "adiabatic_index.h" #include "io_properties.h" #include "kernel_hydro.h" @@ -42,7 +43,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list, list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, UNIT_CONV_LENGTH, parts, h); list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, - UNIT_CONV_ENERGY, parts, entropy); + UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, entropy); list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, parts, id); list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL, @@ -51,6 +52,12 @@ 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) { + + return p->entropy * pow_gamma_minus_one(p->rho) * + hydro_one_over_gamma_minus_one; +} + /** * @brief Specifies which particle fields to write to a dataset * @@ -61,7 +68,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list, void hydro_write_particles(struct part* parts, struct io_props* list, int* num_fields) { - *num_fields = 8; + *num_fields = 9; /* List what we want to write */ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, @@ -80,6 +87,9 @@ 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("InternalEnergy", FLOAT, 1, + UNIT_CONV_ENERGY_PER_UNIT_MASS, + parts, rho, convert_u); } /** @@ -102,4 +112,4 @@ void writeSPHflavour(hid_t h_grpsph) { * * @return 1 if entropy is in 'internal energy', 0 otherwise. */ -int writeEntropyFlag() { return 1; } +int writeEntropyFlag() { return 0; } diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h index 6b52481e9e77ffb7d3c35edd0b2da834ff035a87..1a5b2938354fb8ee3ef936f86ff442c731e0d4c7 100644 --- a/src/hydro/Minimal/hydro_io.h +++ b/src/hydro/Minimal/hydro_io.h @@ -42,7 +42,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list, list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY, UNIT_CONV_LENGTH, parts, h); list[4] = io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY, - UNIT_CONV_ENERGY, parts, u); + UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY, UNIT_CONV_NO_UNITS, parts, id); list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL, @@ -72,7 +72,7 @@ void hydro_write_particles(struct part* parts, struct io_props* list, io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass); list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH, parts, h); - list[4] = io_make_output_field("Entropy", FLOAT, 1, + list[4] = io_make_output_field("InternalEnergy", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u); list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, parts, id); diff --git a/src/io_properties.h b/src/io_properties.h index 5e19420d006ee8043ac418045b6889d885fbd66a..2d4a65163f1ef2bcbaceadee73596d8c0345171a 100644 --- a/src/io_properties.h +++ b/src/io_properties.h @@ -54,6 +54,16 @@ struct io_props { /* The size of the particles */ size_t partSize; + + /* The particle arrays */ + struct part* parts; + struct gpart* gparts; + + /* Conversion function for part */ + float (*convert_part)(struct engine*, struct part*); + + /* Conversion function for gpart */ + float (*convert_gpart)(struct engine*, struct gpart*); }; /** @@ -68,7 +78,7 @@ struct io_props { * * @param name Name of the field to read * @param type The type of the data - * @param dimension Dataset dimension (!D, 3D, ...) + * @param dimension Dataset dimension (1D, 3D, ...) * @param importance Is this dataset compulsory ? * @param units The units of the dataset * @param field Pointer to the field of the first particle @@ -89,6 +99,10 @@ struct io_props io_make_input_field_(char name[FIELD_BUFFER_SIZE], r.units = units; r.field = field; r.partSize = partSize; + r.parts = NULL; + r.gparts = NULL; + r.convert_part = NULL; + r.convert_gpart = NULL; return r; } @@ -105,7 +119,7 @@ struct io_props io_make_input_field_(char name[FIELD_BUFFER_SIZE], * * @param name Name of the field to read * @param type The type of the data - * @param dimension Dataset dimension (!D, 3D, ...) + * @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 @@ -124,6 +138,98 @@ struct io_props io_make_output_field_(char name[FIELD_BUFFER_SIZE], r.units = units; r.field = field; r.partSize = partSize; + r.parts = NULL; + r.gparts = NULL; + r.convert_part = NULL; + r.convert_gpart = NULL; + + return r; +} + +/** + * @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) + +/** + * @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 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_( + char name[FIELD_BUFFER_SIZE], enum DATA_TYPE type, int dimension, + enum UnitConversionFactor units, char* field, size_t partSize, + struct part* parts, float (*functionPtr)(struct engine*, struct part*)) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = 0; + r.units = units; + r.field = field; + r.partSize = partSize; + r.parts = parts; + r.gparts = NULL; + r.convert_part = functionPtr; + r.convert_gpart = NULL; + + return r; +} + +/** + * @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) + +/** + * @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 field Pointer to the field of the first particle + * @param partSize 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_( + char name[FIELD_BUFFER_SIZE], enum DATA_TYPE type, int dimension, + enum UnitConversionFactor units, char* field, size_t partSize, + struct gpart* gparts, float (*functionPtr)(struct engine*, struct gpart*)) { + + struct io_props r; + strcpy(r.name, name); + r.type = type; + r.dimension = dimension; + r.importance = 0; + r.units = units; + r.field = field; + r.partSize = partSize; + r.parts = NULL; + r.gparts = gparts; + r.convert_part = NULL; + r.convert_gpart = functionPtr; return r; } diff --git a/src/parallel_io.c b/src/parallel_io.c index 5e358c34f70065255e1eab61295b21487f07f77d..c6bbcdfd3982ee5663ff54246feb58bebdbb0dad 100644 --- a/src/parallel_io.c +++ b/src/parallel_io.c @@ -181,6 +181,7 @@ void readArray(hid_t grp, const struct io_props prop, size_t N, /** * @brief Writes a data array in given HDF5 group. * + * @param e The #engine we are writing from. * @param grp The group in which to write. * @param fileName The name of the file in which the data is written * @param xmfFile The FILE used to write the XMF description @@ -194,7 +195,7 @@ void readArray(hid_t grp, const struct io_props prop, size_t N, * the part array will be written once the structures have been stabilized. * */ -void writeArray(hid_t grp, char* fileName, FILE* xmfFile, +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, const struct UnitSystem* internal_units, @@ -211,9 +212,25 @@ void writeArray(hid_t grp, char* fileName, FILE* xmfFile, if (temp == NULL) error("Unable to allocate memory for temporary buffer"); /* Copy particle data to temporary buffer */ - char* temp_c = temp; - for (size_t i = 0; i < N; ++i) - memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize); + if (props.convert_part == NULL && + props.convert_gpart == NULL) { /* 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)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_part(e, &props.parts[i]); + + } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_gpart(e, &props.gparts[i]); + } /* Unit conversion if necessary */ const double factor = @@ -798,7 +815,7 @@ void write_output_parallel(struct engine* e, const char* baseName, /* Write everything */ for (int i = 0; i < num_fields; ++i) - writeArray(h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, + writeArray(e, h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, N_total[ptype], mpi_rank, offset[ptype], internal_units, snapshot_units); diff --git a/src/runner.c b/src/runner.c index ae8a33f8b6f62ad0e89e78bc04f60c965c88cd93..4d0b595b8241a8a6887f7da9e657f892c309103f 100644 --- a/src/runner.c +++ b/src/runner.c @@ -408,7 +408,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) { if (gp->ti_end <= ti_current) { /* Get ready for a density calculation */ - gravity_init_part(gp); + gravity_init_gpart(gp); } } } diff --git a/src/serial_io.c b/src/serial_io.c index 0134fa3e52cd0df7341cbc88f11f4182b1c02949..91231e70c603fc882f38e73f3d545828613629fc 100644 --- a/src/serial_io.c +++ b/src/serial_io.c @@ -261,6 +261,7 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, /** * @brief Writes a data array in given HDF5 group. * + * @param e The #engine we are writing from. * @param grp The group in which to write. * @param fileName The name of the file in which the data is written * @param xmfFile The FILE used to write the XMF description @@ -276,7 +277,7 @@ void prepareArray(hid_t grp, char* fileName, FILE* xmfFile, * @param us The UnitSystem currently in use * @param convFactor The UnitConversionFactor for this arrayo */ -void writeArray(hid_t grp, char* fileName, FILE* xmfFile, +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, const struct UnitSystem* internal_units, @@ -298,9 +299,25 @@ void writeArray(hid_t grp, char* fileName, FILE* xmfFile, if (temp == NULL) error("Unable to allocate memory for temporary buffer"); /* Copy particle data to temporary buffer */ - char* temp_c = temp; - for (size_t i = 0; i < N; ++i) - memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize); + if (props.convert_part == NULL && + props.convert_gpart == NULL) { /* 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)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_part(e, &props.parts[i]); + + } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_gpart(e, &props.gparts[i]); + } /* Unit conversion if necessary */ const double factor = @@ -885,7 +902,7 @@ void write_output_serial(struct engine* e, const char* baseName, /* Write everything */ for (int i = 0; i < num_fields; ++i) - writeArray(h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, + writeArray(e, h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, N_total[ptype], mpi_rank, offset[ptype], internal_units, snapshot_units); diff --git a/src/single_io.c b/src/single_io.c index 3f634737194bf8bfdc780e662262b228ebc5c737..f88f3da5c352fc30f1a34de453e9b9abfacd3ad7 100644 --- a/src/single_io.c +++ b/src/single_io.c @@ -152,6 +152,7 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N, /** * @brief Writes a data array in given HDF5 group. * + * @param e The #engine we are writing from. * @param grp The group in which to write. * @param fileName The name of the file in which the data is written * @param xmfFile The FILE used to write the XMF description @@ -165,7 +166,7 @@ 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(hid_t grp, char* fileName, FILE* xmfFile, +void writeArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile, char* partTypeGroupName, const struct io_props props, size_t N, const struct UnitSystem* internal_units, const struct UnitSystem* snapshot_units) { @@ -181,9 +182,25 @@ void writeArray(hid_t grp, char* fileName, FILE* xmfFile, if (temp == NULL) error("Unable to allocate memory for temporary buffer"); /* Copy particle data to temporary buffer */ - char* temp_c = temp; - for (size_t i = 0; i < N; ++i) - memcpy(&temp_c[i * copySize], props.field + i * props.partSize, copySize); + if (props.convert_part == NULL && + props.convert_gpart == NULL) { /* 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)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_part(e, &props.parts[i]); + + } else if (props.convert_gpart != NULL) { /* conversion (for gparts)*/ + + float* temp_f = temp; + for (size_t i = 0; i < N; ++i) + temp_f[i] = props.convert_gpart(e, &props.gparts[i]); + } /* Unit conversion if necessary */ const double factor = @@ -695,7 +712,7 @@ void write_output_single(struct engine* e, const char* baseName, /* Write everything */ for (int i = 0; i < num_fields; ++i) - writeArray(h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, + writeArray(e, h_grp, fileName, xmfFile, partTypeGroupName, list[i], N, internal_units, snapshot_units); /* Free temporary array */ diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c index bca40fedb749e5d01cd7ef8a1fd2c77a6bb93b2a..675d088acfea7bc9414d18995de3101115b5187e 100644 --- a/tests/testSPHStep.c +++ b/tests/testSPHStep.c @@ -96,10 +96,16 @@ void runner_dopair2_force(struct runner *r, struct cell *ci, struct cell *cj); /* Run a full time step integration for one cell */ int main() { +#ifndef DEFAULT_SPH return 0; +#else int i, j, k, offset[3]; struct part *p; + struct hydro_props hp; + hp.target_neighbours = 48.; + hp.delta_neighbours = 1.; + hp.max_smoothing_iterations = 30; int N = 10; float dim = 1.; @@ -144,6 +150,7 @@ int main() { struct engine e; e.s = &space; + e.hydro_properties = &hp; struct runner r; r.e = &e; @@ -197,4 +204,5 @@ int main() { free(ci->xparts); return 0; +#endif }