Commit 7ebaf78e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

The system of units is now carried over all the way to the writing routine

parent 6d9f54ae
......@@ -332,8 +332,8 @@ int main(int argc, char *argv[]) {
myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, dry_run);
#endif
#else
read_ic_single(ICfileName, dim, &parts, &gparts, &Ngas, &Ngpart, &periodic,
dry_run);
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &Ngas, &Ngpart,
&periodic, dry_run);
#endif
if (myrank == 0) {
clocks_gettime(&toc);
......@@ -421,7 +421,7 @@ int main(int argc, char *argv[]) {
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &prog_const, &hydro_properties,
engine_policies, talking, &us, &prog_const, &hydro_properties,
&potential);
if (myrank == 0) {
clocks_gettime(&toc);
......
......@@ -276,12 +276,14 @@ void writeAttribute_s(hid_t grp, const char* name, const char* str) {
/**
* @brief Writes the current Unit System
* @param h_file The (opened) HDF5 file in which to write
* @param us The UnitSystem used in the run
* @param us The UnitSystem to dump
* @param groupName The name of the HDF5 group to write to
*/
void writeUnitSystem(hid_t h_file, struct UnitSystem* us) {
hid_t h_grpunit = 0;
void writeUnitSystem(hid_t h_file, const struct UnitSystem* us,
const char* groupName) {
h_grpunit = H5Gcreate1(h_file, "/Units", 0);
hid_t h_grpunit = 0;
h_grpunit = H5Gcreate1(h_file, groupName, 0);
if (h_grpunit < 0) error("Error while creating Unit System group");
writeAttribute_d(h_grpunit, "Unit mass in cgs (U_M)",
......
......@@ -103,7 +103,8 @@ void writeXMFline(FILE* xmfFile, char* fileName, char* partTypeGroupName,
char* name, size_t N, int dim, enum DATA_TYPE type);
void writeCodeDescription(hid_t h_file);
void writeUnitSystem(hid_t h_file, struct UnitSystem* us);
void writeUnitSystem(hid_t h_grp, const struct UnitSystem* us,
const char* groupName);
#endif /* defined HDF5 */
......
......@@ -2565,7 +2565,8 @@ void engine_dump_snapshot(struct engine *e) {
e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
#endif
#else
write_output_single(e, e->snapshotBaseName, e->snapshotUnits);
write_output_single(e, e->snapshotBaseName, e->internalUnits,
e->snapshotUnits);
#endif
clocks_gettime(&time2);
......@@ -2640,6 +2641,7 @@ void engine_unpin() {
* @param with_aff use processor affinity, if supported.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
* @param internal_units The system of units used internally.
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
......@@ -2647,6 +2649,7 @@ void engine_unpin() {
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
const struct UnitSystem *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential) {
......@@ -2676,6 +2679,7 @@ void engine_init(struct engine *e, struct space *s,
e->timeStep = 0.;
e->timeBase = 0.;
e->timeBase_inv = 0.;
e->internalUnits = internal_units;
e->timeFirstSnapshot =
parser_get_param_double(params, "Snapshots:time_first");
e->deltaTimeSnapshot =
......
......@@ -137,6 +137,9 @@ struct engine {
/* Number of particles updated */
size_t updates, g_updates;
/* The internal system of units */
const struct UnitSystem *internalUnits;
/* Snapshot information */
double timeFirstSnapshot;
double deltaTimeSnapshot;
......@@ -207,6 +210,7 @@ void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int policy, int verbose,
const struct UnitSystem *internal_units,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential);
......
......@@ -57,24 +57,27 @@ __attribute__((always_inline)) INLINE static void darkmatter_read_particles(
* @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
* @param internal_units The #UnitSystem used internally
* @param snapshot_units The #UnitSystem used in the snapshots
*
*/
__attribute__((always_inline)) INLINE static void darkmatter_write_particles(
hid_t h_grp, char* fileName, char* partTypeGroupName, FILE* xmfFile,
int Ndm, long long Ndm_total, int mpi_rank, long long offset,
struct gpart* gparts, struct UnitSystem* us) {
struct gpart* gparts, const struct UnitSystem* internal_units,
const struct UnitSystem* snapshot_units) {
/* Write arrays */
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Coordinates", DOUBLE,
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, x, us,
UNIT_CONV_LENGTH);
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, x, internal_units,
snapshot_units, UNIT_CONV_LENGTH);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Masses", FLOAT, Ndm,
1, gparts, Ndm_total, mpi_rank, offset, mass, us, UNIT_CONV_MASS);
1, gparts, Ndm_total, mpi_rank, offset, mass, internal_units,
snapshot_units, UNIT_CONV_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Velocities", FLOAT,
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, v_full, us,
UNIT_CONV_SPEED);
Ndm, 3, gparts, Ndm_total, mpi_rank, offset, v_full,
internal_units, snapshot_units, UNIT_CONV_SPEED);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset, id, us,
UNIT_CONV_NO_UNITS);
ULONGLONG, Ndm, 1, gparts, Ndm_total, mpi_rank, offset, id,
internal_units, snapshot_units, UNIT_CONV_NO_UNITS);
}
......@@ -65,35 +65,41 @@ __attribute__((always_inline)) INLINE static void hydro_read_particles(
* @param offset The offset of the particles for this MPI rank (only used in MPI
*mode)
* @param parts The particle array
* @param us The unit system to use
* @param internal_units The #UnitSystem used internally
* @param snapshot_units The #UnitSystem used in the snapshots
*
*/
__attribute__((always_inline)) INLINE static void hydro_write_particles(
hid_t h_grp, char* fileName, char* partTypeGroupName, FILE* xmfFile, int N,
long long N_total, int mpi_rank, long long offset, struct part* parts,
struct UnitSystem* us) {
const struct UnitSystem* internal_units,
const struct UnitSystem* snapshot_units) {
/* Write arrays */
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Coordinates", DOUBLE,
N, 3, parts, N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
N, 3, parts, N_total, mpi_rank, offset, x, internal_units,
snapshot_units, UNIT_CONV_LENGTH);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Velocities", FLOAT,
N, 3, parts, N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
N, 3, parts, N_total, mpi_rank, offset, v, internal_units,
snapshot_units, UNIT_CONV_SPEED);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Masses", FLOAT, N, 1,
parts, N_total, mpi_rank, offset, mass, us, UNIT_CONV_MASS);
parts, N_total, mpi_rank, offset, mass, internal_units,
snapshot_units, UNIT_CONV_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "SmoothingLength",
FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, us,
UNIT_CONV_LENGTH);
FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, internal_units,
snapshot_units, UNIT_CONV_LENGTH);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy", FLOAT, N,
1, parts, N_total, mpi_rank, offset, entropy, us,
UNIT_CONV_ENTROPY_PER_UNIT_MASS);
1, parts, N_total, mpi_rank, offset, entropy, internal_units,
snapshot_units, UNIT_CONV_ENTROPY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id, us,
UNIT_CONV_NO_UNITS);
ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id,
internal_units, snapshot_units, UNIT_CONV_NO_UNITS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Acceleration", FLOAT,
N, 3, parts, N_total, mpi_rank, offset, a_hydro, us,
UNIT_CONV_ACCELERATION);
N, 3, parts, N_total, mpi_rank, offset, a_hydro, internal_units,
snapshot_units, UNIT_CONV_ACCELERATION);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Density", FLOAT, N,
1, parts, N_total, mpi_rank, offset, rho, us, UNIT_CONV_DENSITY);
1, parts, N_total, mpi_rank, offset, rho, internal_units,
snapshot_units, UNIT_CONV_DENSITY);
}
/**
......
......@@ -43,6 +43,20 @@
#include "part.h"
#include "units.h"
void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
int* periodic, int mpi_rank, int mpi_size, MPI_Comm comm,
MPI_Info info, int dry_run) {}
void write_output_parallel(struct engine* e, const char* baseName,
struct UnitSystem* us, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info) {}
#endif
#if 0
/**
* @brief Reads a data array from a given HDF5 group.
*
......
......@@ -149,7 +149,8 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
* @param part_c A (char*) pointer on the first occurrence of the field of
*interest in the parts array.
* @param partSize The size in bytes of the particle structure.
* @param us The UnitSystem currently in use
* @param internal_units The #UnitSystem used internally
* @param snapshot_units The #UnitSystem used in the snapshots
* @param convFactor The UnitConversionFactor for this array
*
* @todo A better version using HDF5 hyper-slabs to write the file directly from
......@@ -159,7 +160,8 @@ void readArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
char* partTypeGroupName, char* name, enum DATA_TYPE type,
int N, int dim, char* part_c, size_t partSize,
struct UnitSystem* us,
const struct UnitSystem* internal_units,
const struct UnitSystem* snapshot_units,
enum UnitConversionFactor convFactor) {
hid_t h_data = 0, h_err = 0, h_space = 0, h_prop = 0;
void* temp = 0;
......@@ -244,11 +246,13 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
writeXMFline(xmfFile, fileName, partTypeGroupName, name, N, dim, type);
/* Write unit conversion factors for this data set */
units_conversion_string(buffer, us, convFactor);
units_conversion_string(buffer, snapshot_units, convFactor);
writeAttribute_d(h_data, "CGS conversion factor",
units_conversion_factor(us, convFactor));
writeAttribute_f(h_data, "h-scale exponent", units_h_factor(us, convFactor));
writeAttribute_f(h_data, "a-scale exponent", units_a_factor(us, convFactor));
units_conversion_factor(snapshot_units, convFactor));
writeAttribute_f(h_data, "h-scale exponent",
units_h_factor(snapshot_units, convFactor));
writeAttribute_f(h_data, "a-scale exponent",
units_a_factor(snapshot_units, convFactor));
writeAttribute_s(h_data, "Conversion factor", buffer);
/* Free and close everything */
......@@ -296,16 +300,17 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
* @param mpi_rank Unused parameter in non-MPI mode
* @param offset Unused parameter in non-MPI mode
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
* @param internal_units The #UnitSystem used internally
* @param snapshot_units The #UnitSystem used in the snapshots
* @param convFactor The UnitConversionFactor for this array
*
*/
#define writeArray(grp, fileName, xmfFile, partTypeGroupName, name, type, N, \
dim, part, N_total, mpi_rank, offset, field, us, \
convFactor) \
dim, part, N_total, mpi_rank, offset, field, \
internal_units, snapshot_units, convFactor) \
writeArrayBackEnd(grp, fileName, xmfFile, partTypeGroupName, name, type, N, \
dim, (char*)(&(part[0]).field), sizeof(part[0]), us, \
convFactor)
dim, (char*)(&(part[0]).field), sizeof(part[0]), \
internal_units, snapshot_units, convFactor)
/* Import the right hydro definition */
#include "hydro_io.h"
......@@ -316,6 +321,7 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
* @brief Reads an HDF5 initial condition file (GADGET-3 type)
*
* @param fileName The file to read.
* @param internal_units The system units used internally
* @param dim (output) The dimension of the volume.
* @param parts (output) Array of Gas particles.
* @param gparts (output) Array of #gpart particles.
......@@ -332,9 +338,10 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile,
* @todo Read snapshots distributed in more than one file.
*
*/
void read_ic_single(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ngparts,
int* periodic, int dry_run) {
void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
size_t* Ngas, size_t* Ngparts, int* periodic, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
double boxSize[3] = {0.0, -1.0, -1.0};
......@@ -460,7 +467,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
*
* @param e The engine containing all the system.
* @param baseName The common part of the snapshot file name.
* @param us The UnitSystem used for the conversion of units in the output.
* @param internal_units The #UnitSystem used internally
* @param snapshot_units The #UnitSystem used in the snapshots
*
* Creates an HDF5 output file and writes the particles contained
* in the engine. If such a file already exists, it is erased and replaced
......@@ -471,7 +479,8 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts,
*
*/
void write_output_single(struct engine* e, const char* baseName,
struct UnitSystem* us) {
const struct UnitSystem* internal_units,
const struct UnitSystem* snapshot_units) {
hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
......@@ -572,8 +581,11 @@ void write_output_single(struct engine* e, const char* baseName,
parser_write_params_to_hdf5(e->parameter_file, h_grp);
H5Gclose(h_grp);
/* Print the system of Units */
writeUnitSystem(h_file, us);
/* Print the system of Units used in the spashot */
writeUnitSystem(h_file, snapshot_units, "Units");
/* Print the system of Units used internally */
writeUnitSystem(h_file, internal_units, "InternalUnits");
/* Loop over all particle types */
for (int ptype = 0; ptype < NUM_PARTICLE_TYPES; ptype++) {
......@@ -601,7 +613,8 @@ void write_output_single(struct engine* e, const char* baseName,
case GAS:
hydro_write_particles(h_grp, fileName, partTypeGroupName, xmfFile, Ngas,
Ngas, 0, 0, parts, us);
Ngas, 0, 0, parts, internal_units,
snapshot_units);
break;
case DM:
......@@ -616,7 +629,8 @@ void write_output_single(struct engine* e, const char* baseName,
/* Write DM particles */
darkmatter_write_particles(h_grp, fileName, partTypeGroupName, xmfFile,
Ndm, Ndm, 0, 0, dmparts, us);
Ndm, Ndm, 0, 0, dmparts, internal_units,
snapshot_units);
/* Free temporary array */
free(dmparts);
......
......@@ -29,12 +29,13 @@
#if defined(HAVE_HDF5) && !defined(WITH_MPI)
void read_ic_single(char* fileName, double dim[3], struct part** parts,
struct gpart** gparts, size_t* Ngas, size_t* Ndm,
int* periodic, int dry_run);
void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
size_t* Ngas, size_t* Ndm, int* periodic, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
struct UnitSystem* us);
const struct UnitSystem* internal_units,
const struct UnitSystem* snapshot_units);
#endif
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment