diff --git a/src/Makefile.am b/src/Makefile.am
index 5130264a529d4425d8c49ec964820c604bdc01e3..adabaa7095470b9fec781d2ae586760c80b1c6c6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -111,8 +111,9 @@ AM_SOURCES = space.c space_unique_id.c \
engine_marktasks.c engine_drift.c engine_unskip.c engine_collect_end_of_step.c \
engine_redistribute.c engine_fof.c \
queue.c task.c timers.c debug.c scheduler.c proxy.c version.c \
- parallel_io.c common_io.c single_io.c serial_io.c distributed_io.c xmf.c \
- output_options.c line_of_sight.c restart.c parser.c \
+ common_io.c common_io_copy.c common_io_cells.c common_io_fields.c \
+ single_io.c serial_io.c distributed_io.c parallel_io.c \
+ output_options.c line_of_sight.c restart.c parser.c xmf.c \
kernel_hydro.c tools.c map.c part.c partition.c clocks.c \
physical_constants.c units.c potential.c hydro_properties.c \
threadpool.c cooling.c star_formation.c \
diff --git a/src/common_io.c b/src/common_io.c
index ec0d085d914dffa8093b4a67b0b08647097159c5..081d6edac332f934ef0d914c0796783661cd2b36 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -24,35 +24,24 @@
/* This object's header. */
#include "common_io.h"
-/* Pre-inclusion as needed in other headers */
-#include "engine.h"
-
/* Local includes. */
-#include "black_holes_io.h"
-#include "chemistry_io.h"
-#include "const.h"
-#include "cooling_io.h"
+#include "engine.h"
#include "error.h"
-#include "feedback.h"
-#include "fof_io.h"
-#include "gravity_io.h"
-#include "hydro.h"
-#include "hydro_io.h"
-#include "io_properties.h"
#include "kernel_hydro.h"
#include "part.h"
#include "part_type.h"
-#include "rt_io.h"
-#include "sink_io.h"
-#include "star_formation_io.h"
-#include "stars_io.h"
#include "threadpool.h"
#include "tools.h"
-#include "tracers_io.h"
-#include "units.h"
-#include "velociraptor_io.h"
#include "version.h"
+/* I/O functions of each sub-module */
+#include "chemistry_io.h"
+#include "cooling_io.h"
+#include "feedback.h"
+#include "hydro_io.h"
+#include "stars_io.h"
+#include "tracers_io.h"
+
/* Some standard headers. */
#include <math.h>
#include <stddef.h>
@@ -763,457 +752,6 @@ void io_write_engine_policy(hid_t h_file, const struct engine* e) {
H5Gclose(h_grp);
}
-static long long cell_count_non_inhibited_gas(const struct cell* c) {
- const int total_count = c->hydro.count;
- struct part* parts = c->hydro.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((parts[i].time_bin != time_bin_inhibited) &&
- (parts[i].time_bin != time_bin_not_created)) {
- ++count;
- }
- }
- return count;
-}
-
-static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
- const int total_count = c->grav.count;
- struct gpart* gparts = c->grav.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((gparts[i].time_bin != time_bin_inhibited) &&
- (gparts[i].time_bin != time_bin_not_created) &&
- (gparts[i].type == swift_type_dark_matter)) {
- ++count;
- }
- }
- return count;
-}
-
-static long long cell_count_non_inhibited_background_dark_matter(
- const struct cell* c) {
- const int total_count = c->grav.count;
- struct gpart* gparts = c->grav.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((gparts[i].time_bin != time_bin_inhibited) &&
- (gparts[i].time_bin != time_bin_not_created) &&
- (gparts[i].type == swift_type_dark_matter_background)) {
- ++count;
- }
- }
- return count;
-}
-
-static long long cell_count_non_inhibited_stars(const struct cell* c) {
- const int total_count = c->stars.count;
- struct spart* sparts = c->stars.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((sparts[i].time_bin != time_bin_inhibited) &&
- (sparts[i].time_bin != time_bin_not_created)) {
- ++count;
- }
- }
- return count;
-}
-
-static long long cell_count_non_inhibited_black_holes(const struct cell* c) {
- const int total_count = c->black_holes.count;
- struct bpart* bparts = c->black_holes.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((bparts[i].time_bin != time_bin_inhibited) &&
- (bparts[i].time_bin != time_bin_not_created)) {
- ++count;
- }
- }
- return count;
-}
-
-static long long cell_count_non_inhibited_sinks(const struct cell* c) {
- const int total_count = c->sinks.count;
- struct sink* sinks = c->sinks.parts;
- long long count = 0;
- for (int i = 0; i < total_count; ++i) {
- if ((sinks[i].time_bin != time_bin_inhibited) &&
- (sinks[i].time_bin != time_bin_not_created)) {
- ++count;
- }
- }
- return count;
-}
-
-/**
- * @brief Write a single 1D array to a hdf5 group.
- *
- * This creates a simple Nx1 array with a chunk size of 1024x1.
- * The Fletcher-32 filter is applied to the array.
- *
- * @param h_grp The open hdf5 group.
- * @param n The number of elements in the array.
- * @param array The data to write.
- * @param type The type of the data to write.
- * @param name The name of the array.
- * @param array_content The name of the parent group (only used for error
- * messages).
- */
-void io_write_array(hid_t h_grp, const int n, const void* array,
- const enum IO_DATA_TYPE type, const char* name,
- const char* array_content) {
-
- /* Create memory space */
- const hsize_t shape[2] = {(hsize_t)n, 1};
- hid_t h_space = H5Screate(H5S_SIMPLE);
- if (h_space < 0)
- error("Error while creating data space for %s %s", name, array_content);
- hid_t h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
- if (h_err < 0)
- error("Error while changing shape of %s %s data space.", name,
- array_content);
-
- /* Dataset type */
- hid_t h_type = H5Tcopy(io_hdf5_type(type));
-
- const hsize_t chunk[2] = {(1024 > n ? n : 1024), 1};
- hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
- h_err = H5Pset_chunk(h_prop, 1, chunk);
- if (h_err < 0)
- error("Error while setting chunk shapes of %s %s data space.", name,
- array_content);
-
- /* Impose check-sum to verify data corruption */
- h_err = H5Pset_fletcher32(h_prop);
- if (h_err < 0)
- error("Error while setting check-sum filter on %s %s data space.", name,
- array_content);
-
- /* Write */
- hid_t h_data =
- H5Dcreate(h_grp, name, h_type, h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
- if (h_data < 0)
- error("Error while creating dataspace for %s %s.", name, array_content);
- h_err = H5Dwrite(h_data, io_hdf5_type(type), h_space, H5S_ALL, H5P_DEFAULT,
- array);
- if (h_err < 0) error("Error while writing %s %s.", name, array_content);
- H5Tclose(h_type);
- H5Dclose(h_data);
- H5Pclose(h_prop);
- H5Sclose(h_space);
-}
-
-void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
- const struct cell* cells_top, const int nr_cells,
- const double width[3], const int nodeID,
- const int distributed,
- const long long global_counts[swift_type_count],
- const long long global_offsets[swift_type_count],
- const int num_fields[swift_type_count],
- const struct unit_system* internal_units,
- const struct unit_system* snapshot_units) {
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (distributed) {
- if (global_offsets[0] != 0 || global_offsets[1] != 0 ||
- global_offsets[2] != 0 || global_offsets[3] != 0 ||
- global_offsets[4] != 0 || global_offsets[5] != 0)
- error("Global offset non-zero in the distributed io case!");
- }
-#endif
-
- /* Abort if we don't have any cells yet (i.e. haven't constructed the space)
- */
- if (nr_cells == 0) return;
-
- double cell_width[3] = {width[0], width[1], width[2]};
-
- /* Temporary memory for the cell-by-cell information */
- double* centres = NULL;
- centres = (double*)malloc(3 * nr_cells * sizeof(double));
-
- /* Temporary memory for the cell files ID */
- int* files = NULL;
- files = (int*)malloc(nr_cells * sizeof(int));
-
- /* Count of particles in each cell */
- long long *count_part = NULL, *count_gpart = NULL,
- *count_background_gpart = NULL, *count_spart = NULL,
- *count_bpart = NULL, *count_sink = NULL;
- count_part = (long long*)malloc(nr_cells * sizeof(long long));
- count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
- count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
- count_spart = (long long*)malloc(nr_cells * sizeof(long long));
- count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
- count_sink = (long long*)malloc(nr_cells * sizeof(long long));
-
- /* Global offsets of particles in each cell */
- long long *offset_part = NULL, *offset_gpart = NULL,
- *offset_background_gpart = NULL, *offset_spart = NULL,
- *offset_bpart = NULL, *offset_sink = NULL;
- offset_part = (long long*)malloc(nr_cells * sizeof(long long));
- offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
- offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
- offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
- offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
- offset_sink = (long long*)malloc(nr_cells * sizeof(long long));
-
- /* Offsets of the 0^th element */
- offset_part[0] = 0;
- offset_gpart[0] = 0;
- offset_background_gpart[0] = 0;
- offset_spart[0] = 0;
- offset_bpart[0] = 0;
- offset_sink[0] = 0;
-
- /* Collect the cell information of *local* cells */
- long long local_offset_part = 0;
- long long local_offset_gpart = 0;
- long long local_offset_background_gpart = 0;
- long long local_offset_spart = 0;
- long long local_offset_bpart = 0;
- long long local_offset_sink = 0;
- for (int i = 0; i < nr_cells; ++i) {
-
- /* Store in which file this cell will be found */
- if (distributed) {
- files[i] = cells_top[i].nodeID;
- } else {
- files[i] = 0;
- }
-
- /* Is the cell on this node (i.e. we have full information */
- if (cells_top[i].nodeID == nodeID) {
-
- /* Centre of each cell */
- centres[i * 3 + 0] = cells_top[i].loc[0] + cell_width[0] * 0.5;
- centres[i * 3 + 1] = cells_top[i].loc[1] + cell_width[1] * 0.5;
- centres[i * 3 + 2] = cells_top[i].loc[2] + cell_width[2] * 0.5;
-
- /* Finish by box wrapping to match what is done to the particles */
- centres[i * 3 + 0] = box_wrap(centres[i * 3 + 0], 0.0, dim[0]);
- centres[i * 3 + 1] = box_wrap(centres[i * 3 + 1], 0.0, dim[1]);
- centres[i * 3 + 2] = box_wrap(centres[i * 3 + 2], 0.0, dim[2]);
-
- /* Count real particles that will be written */
- count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
- count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
- count_background_gpart[i] =
- cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
- count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
- count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
- count_sink[i] = cell_count_non_inhibited_sinks(&cells_top[i]);
-
- /* Offsets including the global offset of all particles on this MPI rank
- * Note that in the distributed case, the global offsets are 0 such that
- * we actually compute the offset in the file written by this rank. */
- offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
- offset_gpart[i] =
- local_offset_gpart + global_offsets[swift_type_dark_matter];
- offset_background_gpart[i] =
- local_offset_background_gpart +
- global_offsets[swift_type_dark_matter_background];
- offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
- offset_bpart[i] =
- local_offset_bpart + global_offsets[swift_type_black_hole];
- offset_sink[i] = local_offset_sink + global_offsets[swift_type_sink];
-
- local_offset_part += count_part[i];
- local_offset_gpart += count_gpart[i];
- local_offset_background_gpart += count_background_gpart[i];
- local_offset_spart += count_spart[i];
- local_offset_bpart += count_bpart[i];
- local_offset_sink += count_sink[i];
-
- } else {
-
- /* Just zero everything for the foregin cells */
-
- centres[i * 3 + 0] = 0.;
- centres[i * 3 + 1] = 0.;
- centres[i * 3 + 2] = 0.;
-
- count_part[i] = 0;
- count_gpart[i] = 0;
- count_background_gpart[i] = 0;
- count_spart[i] = 0;
- count_bpart[i] = 0;
- count_sink[i] = 0;
-
- offset_part[i] = 0;
- offset_gpart[i] = 0;
- offset_background_gpart[i] = 0;
- offset_spart[i] = 0;
- offset_bpart[i] = 0;
- offset_sink[i] = 0;
- }
- }
-
-#ifdef WITH_MPI
- /* Now, reduce all the arrays. Note that we use a bit-wise OR here. This
- is safe as we made sure only local cells have non-zero values. */
- MPI_Allreduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
- MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, count_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, count_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
-
- MPI_Allreduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT,
- MPI_BOR, MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
- MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, offset_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
- MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT,
- MPI_BOR, MPI_COMM_WORLD);
- MPI_Allreduce(MPI_IN_PLACE, offset_bpart, nr_cells, MPI_LONG_LONG_INT,
- MPI_BOR, MPI_COMM_WORLD);
-
- /* For the centres we use a sum as MPI does not like bit-wise operations
- on floating point numbers */
- MPI_Allreduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM,
- MPI_COMM_WORLD);
-#endif
-
- /* When writing a single file, only rank 0 writes the meta-data */
- if ((distributed) || (!distributed && nodeID == 0)) {
-
- /* Unit conversion if necessary */
- const double factor = units_conversion_factor(
- internal_units, snapshot_units, UNIT_CONV_LENGTH);
- if (factor != 1.) {
-
- /* Convert the cell centres */
- for (int i = 0; i < nr_cells; ++i) {
- centres[i * 3 + 0] *= factor;
- centres[i * 3 + 1] *= factor;
- centres[i * 3 + 2] *= factor;
- }
-
- /* Convert the cell widths */
- cell_width[0] *= factor;
- cell_width[1] *= factor;
- cell_width[2] *= factor;
- }
-
- /* Write some meta-information first */
- hid_t h_subgrp =
- H5Gcreate(h_grp, "Meta-data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_subgrp < 0) error("Error while creating meta-data sub-group");
- io_write_attribute(h_subgrp, "nr_cells", INT, &nr_cells, 1);
- io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
- io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
- H5Gclose(h_subgrp);
-
- /* Write the centres to the group */
- hsize_t shape[2] = {(hsize_t)nr_cells, 3};
- hid_t h_space = H5Screate(H5S_SIMPLE);
- if (h_space < 0) error("Error while creating data space for cell centres");
- hid_t h_err = H5Sset_extent_simple(h_space, 2, shape, shape);
- if (h_err < 0)
- error("Error while changing shape of gas offsets data space.");
- hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
- H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_data < 0) error("Error while creating dataspace for gas offsets.");
- h_err = H5Dwrite(h_data, io_hdf5_type(DOUBLE), h_space, H5S_ALL,
- H5P_DEFAULT, centres);
- if (h_err < 0) error("Error while writing centres.");
- H5Dclose(h_data);
- H5Sclose(h_space);
-
- /* Group containing the offsets and counts for each particle type */
- hid_t h_grp_offsets = H5Gcreate(h_grp, "OffsetsInFile", H5P_DEFAULT,
- H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp_offsets < 0) error("Error while creating offsets sub-group");
- hid_t h_grp_files =
- H5Gcreate(h_grp, "Files", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp_files < 0) error("Error while creating filess sub-group");
- hid_t h_grp_counts =
- H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp_counts < 0) error("Error while creating counts sub-group");
-
- if (global_counts[swift_type_gas] > 0 && num_fields[swift_type_gas] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType0", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_part, LONGLONG,
- "PartType0", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_part, LONGLONG, "PartType0",
- "counts");
- }
-
- if (global_counts[swift_type_dark_matter] > 0 &&
- num_fields[swift_type_dark_matter] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType1", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_gpart, LONGLONG,
- "PartType1", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_gpart, LONGLONG, "PartType1",
- "counts");
- }
-
- if (global_counts[swift_type_dark_matter_background] > 0 &&
- num_fields[swift_type_dark_matter_background] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType2", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_background_gpart, LONGLONG,
- "PartType2", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_background_gpart, LONGLONG,
- "PartType2", "counts");
- }
-
- if (global_counts[swift_type_sink] > 0 && num_fields[swift_type_sink] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType3", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_sink, LONGLONG,
- "PartType3", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_sink, LONGLONG, "PartType3",
- "counts");
- }
-
- if (global_counts[swift_type_stars] > 0 &&
- num_fields[swift_type_stars] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType4", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_spart, LONGLONG,
- "PartType4", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_spart, LONGLONG, "PartType4",
- "counts");
- }
-
- if (global_counts[swift_type_black_hole] > 0 &&
- num_fields[swift_type_black_hole] > 0) {
- io_write_array(h_grp_files, nr_cells, files, INT, "PartType5", "files");
- io_write_array(h_grp_offsets, nr_cells, offset_bpart, LONGLONG,
- "PartType5", "offsets");
- io_write_array(h_grp_counts, nr_cells, count_bpart, LONGLONG, "PartType5",
- "counts");
- }
-
- H5Gclose(h_grp_offsets);
- H5Gclose(h_grp_files);
- H5Gclose(h_grp_counts);
- }
-
- /* Free everything we allocated */
- free(centres);
- free(files);
- free(count_part);
- free(count_gpart);
- free(count_background_gpart);
- free(count_spart);
- free(count_bpart);
- free(count_sink);
- free(offset_part);
- free(offset_gpart);
- free(offset_background_gpart);
- free(offset_spart);
- free(offset_bpart);
- free(offset_sink);
-}
-
#endif /* HAVE_HDF5 */
/**
@@ -1250,730 +788,6 @@ size_t io_sizeof_type(enum IO_DATA_TYPE type) {
}
}
-/**
- * @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 = (char*)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 xpart* restrict xparts = props.xparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- float* restrict temp_f = (float*)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, xparts + delta + i,
- &temp_f[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #part into a buffer of ints using a
- * conversion function.
- */
-void io_convert_part_i_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 xpart* restrict xparts = props.xparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- int* restrict temp_i = (int*)temp;
- const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_part_i(e, parts + delta + i, xparts + delta + i,
- &temp_i[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 xpart* restrict xparts = props.xparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- double* restrict temp_d = (double*)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, xparts + delta + i,
- &temp_d[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #part into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_part_l_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 xpart* restrict xparts = props.xparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- long long* restrict temp_l = (long long*)temp;
- const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_part_l(e, parts + delta + i, xparts + delta + i,
- &temp_l[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 = (float*)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 ints using a
- * conversion function.
- */
-void io_convert_gpart_i_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? */
- int* restrict temp_i = (int*)temp;
- const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_gpart_i(e, gparts + delta + i, &temp_i[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 = (double*)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 Mapper function to copy #gpart into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_gpart_l_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? */
- long long* restrict temp_l = (long long*)temp;
- const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_gpart_l(e, gparts + delta + i, &temp_l[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #spart into a buffer of floats using a
- * conversion function.
- */
-void io_convert_spart_f_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct spart* restrict sparts = props.sparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- float* restrict temp_f = (float*)temp;
- const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_spart_f(e, sparts + delta + i, &temp_f[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #spart into a buffer of ints using a
- * conversion function.
- */
-void io_convert_spart_i_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct spart* restrict sparts = props.sparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- int* restrict temp_i = (int*)temp;
- const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_spart_i(e, sparts + delta + i, &temp_i[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #spart into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_spart_d_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct spart* restrict sparts = props.sparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- double* restrict temp_d = (double*)temp;
- const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_spart_d(e, sparts + delta + i, &temp_d[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #spart into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_spart_l_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct spart* restrict sparts = props.sparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- long long* restrict temp_l = (long long*)temp;
- const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_spart_l(e, sparts + delta + i, &temp_l[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #bpart into a buffer of floats using a
- * conversion function.
- */
-void io_convert_bpart_f_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct bpart* restrict bparts = props.bparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- float* restrict temp_f = (float*)temp;
- const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_bpart_f(e, bparts + delta + i, &temp_f[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #bpart into a buffer of ints using a
- * conversion function.
- */
-void io_convert_bpart_i_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct bpart* restrict bparts = props.bparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- int* restrict temp_i = (int*)temp;
- const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_bpart_i(e, bparts + delta + i, &temp_i[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #bpart into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_bpart_d_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct bpart* restrict bparts = props.bparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- double* restrict temp_d = (double*)temp;
- const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_bpart_d(e, bparts + delta + i, &temp_d[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #bpart into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_bpart_l_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct bpart* restrict bparts = props.bparts;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- long long* restrict temp_l = (long long*)temp;
- const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_bpart_l(e, bparts + delta + i, &temp_l[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #sink into a buffer of floats using a
- * conversion function.
- */
-void io_convert_sink_f_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct sink* restrict sinks = props.sinks;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- float* restrict temp_f = (float*)temp;
- const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_sink_f(e, sinks + delta + i, &temp_f[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #sink into a buffer of ints using a
- * conversion function.
- */
-void io_convert_sink_i_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct sink* restrict sinks = props.sinks;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- int* restrict temp_i = (int*)temp;
- const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_sink_i(e, sinks + delta + i, &temp_i[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #sink into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_sink_d_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct sink* restrict sinks = props.sinks;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- double* restrict temp_d = (double*)temp;
- const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_sink_d(e, sinks + delta + i, &temp_d[i * dim]);
-}
-
-/**
- * @brief Mapper function to copy #sink into a buffer of doubles using a
- * conversion function.
- */
-void io_convert_sink_l_mapper(void* restrict temp, int N,
- void* restrict extra_data) {
-
- const struct io_props props = *((const struct io_props*)extra_data);
- const struct sink* restrict sinks = props.sinks;
- const struct engine* e = props.e;
- const size_t dim = props.dimension;
-
- /* How far are we with this chunk? */
- long long* restrict temp_l = (long long*)temp;
- const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
-
- for (int i = 0; i < N; i++)
- props.convert_sink_l(e, sinks + delta + i, &temp_l[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 = (char*)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, threadpool_auto_chunk_size, (void*)&props);
-
- } else { /* Converting particle to data */
-
- if (props.convert_part_f != NULL) {
-
- /* Prepare some parameters */
- float* temp_f = (float*)temp;
- props.start_temp_f = (float*)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,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_part_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_part_i_mapper, temp_i, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_part_d != NULL) {
-
- /* Prepare some parameters */
- double* temp_d = (double*)temp;
- props.start_temp_d = (double*)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,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_part_l != NULL) {
-
- /* Prepare some parameters */
- long long* temp_l = (long long*)temp;
- props.start_temp_l = (long long*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_part_l_mapper, temp_l, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_gpart_f != NULL) {
-
- /* Prepare some parameters */
- float* temp_f = (float*)temp;
- props.start_temp_f = (float*)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,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_gpart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_gpart_i_mapper, temp_i, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_gpart_d != NULL) {
-
- /* Prepare some parameters */
- double* temp_d = (double*)temp;
- props.start_temp_d = (double*)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,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_gpart_l != NULL) {
-
- /* Prepare some parameters */
- long long* temp_l = (long long*)temp;
- props.start_temp_l = (long long*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_gpart_l_mapper, temp_l, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_spart_f != NULL) {
-
- /* Prepare some parameters */
- float* temp_f = (float*)temp;
- props.start_temp_f = (float*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_f_mapper, temp_f, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_spart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_i_mapper, temp_i, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_spart_d != NULL) {
-
- /* Prepare some parameters */
- double* temp_d = (double*)temp;
- props.start_temp_d = (double*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_d_mapper, temp_d, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_spart_l != NULL) {
-
- /* Prepare some parameters */
- long long* temp_l = (long long*)temp;
- props.start_temp_l = (long long*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_spart_l_mapper, temp_l, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_sink_f != NULL) {
-
- /* Prepare some parameters */
- float* temp_f = (float*)temp;
- props.start_temp_f = (float*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_sink_f_mapper, temp_f, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_sink_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_sink_i_mapper, temp_i, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_sink_d != NULL) {
-
- /* Prepare some parameters */
- double* temp_d = (double*)temp;
- props.start_temp_d = (double*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_sink_d_mapper, temp_d, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_sink_l != NULL) {
-
- /* Prepare some parameters */
- long long* temp_l = (long long*)temp;
- props.start_temp_l = (long long*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_sink_l_mapper, temp_l, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_bpart_f != NULL) {
-
- /* Prepare some parameters */
- float* temp_f = (float*)temp;
- props.start_temp_f = (float*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_bpart_f_mapper, temp_f, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_bpart_i != NULL) {
-
- /* Prepare some parameters */
- int* temp_i = (int*)temp;
- props.start_temp_i = (int*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_bpart_i_mapper, temp_i, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_bpart_d != NULL) {
-
- /* Prepare some parameters */
- double* temp_d = (double*)temp;
- props.start_temp_d = (double*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_bpart_d_mapper, temp_d, N, copySize,
- threadpool_auto_chunk_size, (void*)&props);
-
- } else if (props.convert_bpart_l != NULL) {
-
- /* Prepare some parameters */
- long long* temp_l = (long long*)temp;
- props.start_temp_l = (long long*)temp;
- props.e = e;
-
- /* Copy the whole thing into a buffer */
- threadpool_map((struct threadpool*)&e->threadpool,
- io_convert_bpart_l_mapper, temp_l, N, copySize,
- threadpool_auto_chunk_size, (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, (double*)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, (float*)temp,
- IO_BUFFER_ALIGNMENT);
- for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
- }
- }
-}
-
void io_prepare_dm_gparts_mapper(void* restrict data, int Ndm, void* dummy) {
struct gpart* restrict gparts = (struct gpart*)data;
@@ -2536,243 +1350,6 @@ void io_collect_gparts_background_to_write(
count, Ngparts_written);
}
-/**
- * @brief Prepare the output option fields according to the user's choices and
- * verify that they are valid.
- *
- * @param output_options The #output_options for this run
- * @param with_cosmology Ran with cosmology?
- * @param with_fof Are we running with on-the-fly Fof?
- * @param with_stf Are we running with on-the-fly structure finder?
- */
-void io_prepare_output_fields(struct output_options* output_options,
- const int with_cosmology, const int with_fof,
- const int with_stf) {
-
- const int MAX_NUM_PTYPE_FIELDS = 100;
-
- /* Parameter struct for the output options */
- struct swift_params* params = output_options->select_output;
-
- /* Get all possible outputs per particle type */
- int ptype_num_fields_total[swift_type_count] = {0};
- struct io_props field_list[swift_type_count][MAX_NUM_PTYPE_FIELDS];
-
- for (int ptype = 0; ptype < swift_type_count; ptype++)
- ptype_num_fields_total[ptype] = get_ptype_fields(
- ptype, field_list[ptype], with_cosmology, with_fof, with_stf);
-
- /* Check for whether we have a `Default` section */
- int have_default = 0;
-
- /* Loop over each section, i.e. different class of output */
- for (int section_id = 0; section_id < params->sectionCount; section_id++) {
-
- /* Get the name of current (selection) section, without a trailing colon */
- char section_name[FIELD_BUFFER_SIZE];
- strcpy(section_name, params->section[section_id].name);
- section_name[strlen(section_name) - 1] = 0;
-
- /* Is this the `Default` section? */
- if (strcmp(section_name, select_output_header_default_name) == 0)
- have_default = 1;
-
- /* How many fields should each ptype write by default? */
- int ptype_num_fields_to_write[swift_type_count];
-
- /* What is the default writing status for each ptype (on/off)? */
- int ptype_default_write_status[swift_type_count];
-
- /* Initialise section-specific writing counters for each particle type.
- * If default is 'write', then we start from the total to deduct any fields
- * that are switched off. If the default is 'off', we have to start from
- * zero and then count upwards for each field that is switched back on. */
- for (int ptype = 0; ptype < swift_type_count; ptype++) {
-
- /* Internally also verifies that the default level is allowed */
- const enum lossy_compression_schemes compression_level_current_default =
- output_options_get_ptype_default_compression(params, section_name,
- (enum part_type)ptype);
-
- if (compression_level_current_default == compression_do_not_write) {
- ptype_default_write_status[ptype] = 0;
- ptype_num_fields_to_write[ptype] = 0;
- } else {
- ptype_default_write_status[ptype] = 1;
- ptype_num_fields_to_write[ptype] = ptype_num_fields_total[ptype];
- }
-
- } /* ends loop over particle types */
-
- /* Loop over each parameter */
- for (int param_id = 0; param_id < params->paramCount; param_id++) {
-
- /* Full name of the parameter to check */
- const char* param_name = params->data[param_id].name;
-
- /* Check whether the file still contains the old, now inappropriate
- * 'SelectOutput' section */
- if (strstr(param_name, "SelectOutput:") != NULL) {
- error(
- "Output selection files no longer require the use of top level "
- "SelectOutput; see the documentation for changes.");
- }
-
- /* Skip if the parameter belongs to another output class or is a
- * 'Standard' parameter */
- if (strstr(param_name, section_name) == NULL) continue;
- if (strstr(param_name, ":Standard_") != NULL) continue;
-
- /* Get the particle type for current parameter
- * (raises an error if it could not determine it) */
- const int param_ptype = get_param_ptype(param_name);
-
- /* Issue a warning if this parameter does not pertain to any of the
- * known fields from this ptype. */
- int field_id = 0;
- char field_name[PARSER_MAX_LINE_SIZE];
- for (field_id = 0; field_id < ptype_num_fields_total[param_ptype];
- field_id++) {
-
- sprintf(field_name, "%s:%.*s_%s", section_name, FIELD_BUFFER_SIZE,
- field_list[param_ptype][field_id].name,
- part_type_names[param_ptype]);
-
- if (strcmp(param_name, field_name) == 0) break;
- }
-
- int param_is_known = 0; /* Update below if it is a known one */
- if (field_id < ptype_num_fields_total[param_ptype])
- param_is_known = 1;
- else
- message(
- "WARNING: Trying to change behaviour of field '%s' (read from "
- "'%s') that does not exist. This may be because you are not "
- "running with all of the physics that you compiled the code with.",
- param_name, params->fileName);
-
- /* Perform a correctness check on the _value_ of the parameter */
- char param_value[FIELD_BUFFER_SIZE];
- parser_get_param_string(params, param_name, param_value);
-
- int value_id = 0;
- for (value_id = 0; value_id < compression_level_count; value_id++)
- if (strcmp(param_value, lossy_compression_schemes_names[value_id]) == 0)
- break;
-
- if (value_id == compression_level_count)
- error("Choice of output selection parameter %s ('%s') is invalid.",
- param_name, param_value);
-
- /* Adjust number of fields to be written for param_ptype, if this field's
- * status is different from default and it is a known one. */
- if (param_is_known) {
- const int is_on =
- strcmp(param_value,
- lossy_compression_schemes_names[compression_do_not_write]) !=
- 0;
-
- if (is_on && !ptype_default_write_status[param_ptype]) {
- /* Particle should be written even though default is off:
- * increase field count */
- ptype_num_fields_to_write[param_ptype] += 1;
- }
- if (!is_on && ptype_default_write_status[param_ptype]) {
- /* Particle should not be written, even though default is on:
- * decrease field count */
- ptype_num_fields_to_write[param_ptype] -= 1;
- }
- }
- } /* ends loop over parameters */
-
- /* Second loop over ptypes, to write out total number of fields to write */
- for (int ptype = 0; ptype < swift_type_count; ptype++) {
-
-#ifdef SWIFT_DEBUG_CHECKS
- /* Sanity check: is the number of fields to write non-negative? */
- if (ptype_num_fields_to_write[ptype] < 0)
- error(
- "We seem to have subtracted too many fields for particle "
- "type %d in output class %s (total to write is %d)",
- ptype, section_name, ptype_num_fields_to_write[ptype]);
-#endif
- output_options->num_fields_to_write[section_id][ptype] =
- ptype_num_fields_to_write[ptype];
- }
- } /* Ends loop over sections, for different output classes */
-
- /* Add field numbers for (possible) implicit `Default` output class */
- if (!have_default) {
- const int default_id = output_options->select_output->sectionCount;
- for (int ptype = 0; ptype < swift_type_count; ptype++)
- output_options->num_fields_to_write[default_id][ptype] =
- ptype_num_fields_total[ptype];
- }
-}
-
-/**
- * @brief Write the output field parameters file
- *
- * @param filename The file to write.
- * @param with_cosmology Use cosmological name variant?
- */
-void io_write_output_field_parameter(const char* filename, int with_cosmology) {
-
- FILE* file = fopen(filename, "w");
- if (file == NULL) error("Error opening file '%s'", filename);
-
- /* Create a fake unit system for the snapshots */
- struct unit_system snapshot_units;
- units_init_cgs(&snapshot_units);
-
- /* Loop over all particle types */
- fprintf(file, "Default:\n");
- for (int ptype = 0; ptype < swift_type_count; ptype++) {
-
- struct io_props list[100];
- int num_fields = get_ptype_fields(ptype, list, with_cosmology,
- /*with_fof=*/1, /*with_stf=*/1);
-
- if (num_fields == 0) continue;
-
- /* Output a header for that particle type */
- fprintf(file, " # Particle Type %s\n", part_type_names[ptype]);
-
- /* Write all the fields of this particle type */
- for (int i = 0; i < num_fields; ++i) {
-
- char unit_buffer[FIELD_BUFFER_SIZE] = {0};
- units_cgs_conversion_string(unit_buffer, &snapshot_units, list[i].units,
- list[i].scale_factor_exponent);
-
- /* Need to buffer with a maximal size - otherwise we can't read in again
- * because comments are too long */
- char comment_write_buffer[PARSER_MAX_LINE_SIZE / 2];
-
- sprintf(comment_write_buffer, "%.*s", PARSER_MAX_LINE_SIZE / 2 - 1,
- list[i].description);
-
- /* If our string is too long, replace the last few characters (before
- * \0) with ... for 'fancy printing' */
- if (strlen(comment_write_buffer) > PARSER_MAX_LINE_SIZE / 2 - 3) {
- strcpy(&comment_write_buffer[PARSER_MAX_LINE_SIZE / 2 - 4], "...");
- }
-
- fprintf(file, " %s_%s: %s # %s : %s\n", list[i].name,
- part_type_names[ptype], "on", comment_write_buffer, unit_buffer);
- }
-
- fprintf(file, "\n");
- }
-
- fclose(file);
-
- printf(
- "List of valid ouput fields for the particle in snapshots dumped in "
- "'%s'.\n",
- filename);
-}
-
/**
* @brief Create the subdirectory for snapshots if the user demanded one.
*
@@ -2831,114 +1408,6 @@ void io_get_snapshot_filename(char filename[1024], char xmf_filename[1024],
sprintf(xmf_filename, "%s.xmf", basename);
}
}
-
-/**
- * @brief Return the number and names of all output fields of a given ptype.
- *
- * @param ptype The index of the particle type under consideration.
- * @param list An io_props list that will hold the individual fields.
- * @param with_cosmology Use cosmological name variant?
- * @param with_fof Include FoF related fields?
- * @param with_stf Include STF related fields?
- *
- * @return The total number of fields that can be written for the ptype.
- */
-int get_ptype_fields(const int ptype, struct io_props* list,
- const int with_cosmology, const int with_fof,
- const int with_stf) {
-
- int num_fields = 0;
-
- switch (ptype) {
-
- case swift_type_gas:
- hydro_write_particles(NULL, NULL, list, &num_fields);
- num_fields += chemistry_write_particles(NULL, NULL, list + num_fields,
- with_cosmology);
- num_fields +=
- cooling_write_particles(NULL, NULL, list + num_fields, NULL);
- num_fields += tracers_write_particles(NULL, NULL, list + num_fields,
- with_cosmology);
- num_fields +=
- star_formation_write_particles(NULL, NULL, list + num_fields);
- if (with_fof)
- num_fields += fof_write_parts(NULL, NULL, list + num_fields);
- if (with_stf)
- num_fields += velociraptor_write_parts(NULL, NULL, list + num_fields);
- num_fields += rt_write_particles(NULL, list + num_fields);
- break;
-
- case swift_type_dark_matter:
- darkmatter_write_particles(NULL, list, &num_fields);
- if (with_fof) num_fields += fof_write_gparts(NULL, list + num_fields);
- if (with_stf)
- num_fields += velociraptor_write_gparts(NULL, list + num_fields);
- break;
-
- case swift_type_dark_matter_background:
- darkmatter_write_particles(NULL, list, &num_fields);
- if (with_fof) num_fields += fof_write_gparts(NULL, list + num_fields);
- if (with_stf)
- num_fields += velociraptor_write_gparts(NULL, list + num_fields);
- break;
-
- case swift_type_stars:
- stars_write_particles(NULL, list, &num_fields, with_cosmology);
- num_fields += chemistry_write_sparticles(NULL, list + num_fields);
- num_fields +=
- tracers_write_sparticles(NULL, list + num_fields, with_cosmology);
- num_fields += star_formation_write_sparticles(NULL, list + num_fields);
- if (with_fof) num_fields += fof_write_sparts(NULL, list + num_fields);
- if (with_stf)
- num_fields += velociraptor_write_sparts(NULL, list + num_fields);
- num_fields += rt_write_stars(NULL, list + num_fields);
- break;
-
- case swift_type_sink:
- sink_write_particles(NULL, list, &num_fields, with_cosmology);
- break;
-
- case swift_type_black_hole:
- black_holes_write_particles(NULL, list, &num_fields, with_cosmology);
- num_fields += chemistry_write_bparticles(NULL, list + num_fields);
- if (with_fof) num_fields += fof_write_bparts(NULL, list + num_fields);
- if (with_stf)
- num_fields += velociraptor_write_bparts(NULL, list + num_fields);
- break;
-
- default:
- error("Particle Type %d not yet supported. Aborting", ptype);
- }
-
- return num_fields;
-}
-
-/**
- * @brief Return the particle type code of a select_output parameter
- *
- * @param name The name of the parameter under consideration.
- *
- * @return The (integer) particle type of the parameter.
- */
-int get_param_ptype(const char* name) {
-
- const int name_len = strlen(name);
-
- for (int ptype = 0; ptype < swift_type_count; ptype++) {
- const int ptype_name_len = strlen(part_type_names[ptype]);
- if (name_len >= ptype_name_len &&
- strcmp(&name[name_len - ptype_name_len], part_type_names[ptype]) == 0)
- return ptype;
- }
-
- /* If we get here, we could not match the name, so something's gone wrong. */
- error("Could not determine the particle type for parameter '%s'.", name);
-
- /* We can never get here, but the compiler may complain if we don't return
- * an int after promising to do so... */
- return -1;
-}
-
/**
* @brief Set all ParticleIDs for each gpart to 1.
*
@@ -2950,6 +1419,6 @@ int get_param_ptype(const char* name) {
* @param gparts The array of loaded gparts.
* @param Ngparts Number of loaded gparts.
*/
-void set_ids_to_one(struct gpart* gparts, const size_t Ngparts) {
+void io_set_ids_to_one(struct gpart* gparts, const size_t Ngparts) {
for (size_t i = 0; i < Ngparts; i++) gparts[i].id_or_neg_offset = 1;
}
diff --git a/src/common_io.h b/src/common_io.h
index e8fc1fc4f68fbdcae027974593b6d4f124cd1fc7..95f00cf6c3426cb99b6c1d045dd265302cbad468 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -193,10 +193,6 @@ void io_get_snapshot_filename(char filename[1024], char xmf_filename[1024],
const int stf_count, const int snap_count,
const char* subdir, const char* basename);
-int get_ptype_fields(const int ptype, struct io_props* list,
- const int with_cosmology, const int with_fof,
- const int with_stf);
-int get_param_ptype(const char* name);
-void set_ids_to_one(struct gpart* gparts, const size_t Ngparts);
+void io_set_ids_to_one(struct gpart* gparts, const size_t Ngparts);
#endif /* SWIFT_COMMON_IO_H */
diff --git a/src/common_io_cells.c b/src/common_io_cells.c
new file mode 100644
index 0000000000000000000000000000000000000000..31bc6521e80cf5ae4379e3df7f96e474d0e172b0
--- /dev/null
+++ b/src/common_io_cells.c
@@ -0,0 +1,491 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Matthieu Schaller (schaller@strw.leidenuniv.nl)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "common_io.h"
+
+/* Local includes. */
+#include "cell.h"
+#include "timeline.h"
+#include "units.h"
+
+#if defined(HAVE_HDF5)
+
+#include <hdf5.h>
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
+static long long cell_count_non_inhibited_gas(const struct cell* c) {
+ const int total_count = c->hydro.count;
+ struct part* parts = c->hydro.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((parts[i].time_bin != time_bin_inhibited) &&
+ (parts[i].time_bin != time_bin_not_created)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+static long long cell_count_non_inhibited_dark_matter(const struct cell* c) {
+ const int total_count = c->grav.count;
+ struct gpart* gparts = c->grav.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((gparts[i].time_bin != time_bin_inhibited) &&
+ (gparts[i].time_bin != time_bin_not_created) &&
+ (gparts[i].type == swift_type_dark_matter)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+static long long cell_count_non_inhibited_background_dark_matter(
+ const struct cell* c) {
+ const int total_count = c->grav.count;
+ struct gpart* gparts = c->grav.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((gparts[i].time_bin != time_bin_inhibited) &&
+ (gparts[i].time_bin != time_bin_not_created) &&
+ (gparts[i].type == swift_type_dark_matter_background)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+static long long cell_count_non_inhibited_stars(const struct cell* c) {
+ const int total_count = c->stars.count;
+ struct spart* sparts = c->stars.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((sparts[i].time_bin != time_bin_inhibited) &&
+ (sparts[i].time_bin != time_bin_not_created)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+static long long cell_count_non_inhibited_black_holes(const struct cell* c) {
+ const int total_count = c->black_holes.count;
+ struct bpart* bparts = c->black_holes.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((bparts[i].time_bin != time_bin_inhibited) &&
+ (bparts[i].time_bin != time_bin_not_created)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+static long long cell_count_non_inhibited_sinks(const struct cell* c) {
+ const int total_count = c->sinks.count;
+ struct sink* sinks = c->sinks.parts;
+ long long count = 0;
+ for (int i = 0; i < total_count; ++i) {
+ if ((sinks[i].time_bin != time_bin_inhibited) &&
+ (sinks[i].time_bin != time_bin_not_created)) {
+ ++count;
+ }
+ }
+ return count;
+}
+
+/**
+ * @brief Write a single 1D array to a hdf5 group.
+ *
+ * This creates a simple Nx1 array with a chunk size of 1024x1.
+ * The Fletcher-32 filter is applied to the array.
+ *
+ * @param h_grp The open hdf5 group.
+ * @param n The number of elements in the array.
+ * @param array The data to write.
+ * @param type The type of the data to write.
+ * @param name The name of the array.
+ * @param array_content The name of the parent group (only used for error
+ * messages).
+ */
+void io_write_array(hid_t h_grp, const int n, const void* array,
+ const enum IO_DATA_TYPE type, const char* name,
+ const char* array_content) {
+
+ /* Create memory space */
+ const hsize_t shape[2] = {(hsize_t)n, 1};
+ hid_t h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for %s %s", name, array_content);
+ hid_t h_err = H5Sset_extent_simple(h_space, 1, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of %s %s data space.", name,
+ array_content);
+
+ /* Dataset type */
+ hid_t h_type = H5Tcopy(io_hdf5_type(type));
+
+ const hsize_t chunk[2] = {(1024 > n ? n : 1024), 1};
+ hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
+ h_err = H5Pset_chunk(h_prop, 1, chunk);
+ if (h_err < 0)
+ error("Error while setting chunk shapes of %s %s data space.", name,
+ array_content);
+
+ /* Impose check-sum to verify data corruption */
+ h_err = H5Pset_fletcher32(h_prop);
+ if (h_err < 0)
+ error("Error while setting check-sum filter on %s %s data space.", name,
+ array_content);
+
+ /* Write */
+ hid_t h_data =
+ H5Dcreate(h_grp, name, h_type, h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
+ if (h_data < 0)
+ error("Error while creating dataspace for %s %s.", name, array_content);
+ h_err = H5Dwrite(h_data, io_hdf5_type(type), h_space, H5S_ALL, H5P_DEFAULT,
+ array);
+ if (h_err < 0) error("Error while writing %s %s.", name, array_content);
+ H5Tclose(h_type);
+ H5Dclose(h_data);
+ H5Pclose(h_prop);
+ H5Sclose(h_space);
+}
+
+void io_write_cell_offsets(hid_t h_grp, const int cdim[3], const double dim[3],
+ const struct cell* cells_top, const int nr_cells,
+ const double width[3], const int nodeID,
+ const int distributed,
+ const long long global_counts[swift_type_count],
+ const long long global_offsets[swift_type_count],
+ const int num_fields[swift_type_count],
+ const struct unit_system* internal_units,
+ const struct unit_system* snapshot_units) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (distributed) {
+ if (global_offsets[0] != 0 || global_offsets[1] != 0 ||
+ global_offsets[2] != 0 || global_offsets[3] != 0 ||
+ global_offsets[4] != 0 || global_offsets[5] != 0)
+ error("Global offset non-zero in the distributed io case!");
+ }
+#endif
+
+ /* Abort if we don't have any cells yet (i.e. haven't constructed the space)
+ */
+ if (nr_cells == 0) return;
+
+ double cell_width[3] = {width[0], width[1], width[2]};
+
+ /* Temporary memory for the cell-by-cell information */
+ double* centres = NULL;
+ centres = (double*)malloc(3 * nr_cells * sizeof(double));
+
+ /* Temporary memory for the cell files ID */
+ int* files = NULL;
+ files = (int*)malloc(nr_cells * sizeof(int));
+
+ /* Count of particles in each cell */
+ long long *count_part = NULL, *count_gpart = NULL,
+ *count_background_gpart = NULL, *count_spart = NULL,
+ *count_bpart = NULL, *count_sink = NULL;
+ count_part = (long long*)malloc(nr_cells * sizeof(long long));
+ count_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ count_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ count_spart = (long long*)malloc(nr_cells * sizeof(long long));
+ count_bpart = (long long*)malloc(nr_cells * sizeof(long long));
+ count_sink = (long long*)malloc(nr_cells * sizeof(long long));
+
+ /* Global offsets of particles in each cell */
+ long long *offset_part = NULL, *offset_gpart = NULL,
+ *offset_background_gpart = NULL, *offset_spart = NULL,
+ *offset_bpart = NULL, *offset_sink = NULL;
+ offset_part = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_background_gpart = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_spart = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_bpart = (long long*)malloc(nr_cells * sizeof(long long));
+ offset_sink = (long long*)malloc(nr_cells * sizeof(long long));
+
+ /* Offsets of the 0^th element */
+ offset_part[0] = 0;
+ offset_gpart[0] = 0;
+ offset_background_gpart[0] = 0;
+ offset_spart[0] = 0;
+ offset_bpart[0] = 0;
+ offset_sink[0] = 0;
+
+ /* Collect the cell information of *local* cells */
+ long long local_offset_part = 0;
+ long long local_offset_gpart = 0;
+ long long local_offset_background_gpart = 0;
+ long long local_offset_spart = 0;
+ long long local_offset_bpart = 0;
+ long long local_offset_sink = 0;
+ for (int i = 0; i < nr_cells; ++i) {
+
+ /* Store in which file this cell will be found */
+ if (distributed) {
+ files[i] = cells_top[i].nodeID;
+ } else {
+ files[i] = 0;
+ }
+
+ /* Is the cell on this node (i.e. we have full information */
+ if (cells_top[i].nodeID == nodeID) {
+
+ /* Centre of each cell */
+ centres[i * 3 + 0] = cells_top[i].loc[0] + cell_width[0] * 0.5;
+ centres[i * 3 + 1] = cells_top[i].loc[1] + cell_width[1] * 0.5;
+ centres[i * 3 + 2] = cells_top[i].loc[2] + cell_width[2] * 0.5;
+
+ /* Finish by box wrapping to match what is done to the particles */
+ centres[i * 3 + 0] = box_wrap(centres[i * 3 + 0], 0.0, dim[0]);
+ centres[i * 3 + 1] = box_wrap(centres[i * 3 + 1], 0.0, dim[1]);
+ centres[i * 3 + 2] = box_wrap(centres[i * 3 + 2], 0.0, dim[2]);
+
+ /* Count real particles that will be written */
+ count_part[i] = cell_count_non_inhibited_gas(&cells_top[i]);
+ count_gpart[i] = cell_count_non_inhibited_dark_matter(&cells_top[i]);
+ count_background_gpart[i] =
+ cell_count_non_inhibited_background_dark_matter(&cells_top[i]);
+ count_spart[i] = cell_count_non_inhibited_stars(&cells_top[i]);
+ count_bpart[i] = cell_count_non_inhibited_black_holes(&cells_top[i]);
+ count_sink[i] = cell_count_non_inhibited_sinks(&cells_top[i]);
+
+ /* Offsets including the global offset of all particles on this MPI rank
+ * Note that in the distributed case, the global offsets are 0 such that
+ * we actually compute the offset in the file written by this rank. */
+ offset_part[i] = local_offset_part + global_offsets[swift_type_gas];
+ offset_gpart[i] =
+ local_offset_gpart + global_offsets[swift_type_dark_matter];
+ offset_background_gpart[i] =
+ local_offset_background_gpart +
+ global_offsets[swift_type_dark_matter_background];
+ offset_spart[i] = local_offset_spart + global_offsets[swift_type_stars];
+ offset_bpart[i] =
+ local_offset_bpart + global_offsets[swift_type_black_hole];
+ offset_sink[i] = local_offset_sink + global_offsets[swift_type_sink];
+
+ local_offset_part += count_part[i];
+ local_offset_gpart += count_gpart[i];
+ local_offset_background_gpart += count_background_gpart[i];
+ local_offset_spart += count_spart[i];
+ local_offset_bpart += count_bpart[i];
+ local_offset_sink += count_sink[i];
+
+ } else {
+
+ /* Just zero everything for the foregin cells */
+
+ centres[i * 3 + 0] = 0.;
+ centres[i * 3 + 1] = 0.;
+ centres[i * 3 + 2] = 0.;
+
+ count_part[i] = 0;
+ count_gpart[i] = 0;
+ count_background_gpart[i] = 0;
+ count_spart[i] = 0;
+ count_bpart[i] = 0;
+ count_sink[i] = 0;
+
+ offset_part[i] = 0;
+ offset_gpart[i] = 0;
+ offset_background_gpart[i] = 0;
+ offset_spart[i] = 0;
+ offset_bpart[i] = 0;
+ offset_sink[i] = 0;
+ }
+ }
+
+#ifdef WITH_MPI
+ /* Now, reduce all the arrays. Note that we use a bit-wise OR here. This
+ is safe as we made sure only local cells have non-zero values. */
+ MPI_Allreduce(MPI_IN_PLACE, count_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, count_gpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, count_background_gpart, nr_cells,
+ MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, count_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, count_spart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, count_bpart, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+
+ MPI_Allreduce(MPI_IN_PLACE, offset_part, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, offset_gpart, nr_cells, MPI_LONG_LONG_INT,
+ MPI_BOR, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, offset_background_gpart, nr_cells,
+ MPI_LONG_LONG_INT, MPI_BOR, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, offset_sink, nr_cells, MPI_LONG_LONG_INT, MPI_BOR,
+ MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, offset_spart, nr_cells, MPI_LONG_LONG_INT,
+ MPI_BOR, MPI_COMM_WORLD);
+ MPI_Allreduce(MPI_IN_PLACE, offset_bpart, nr_cells, MPI_LONG_LONG_INT,
+ MPI_BOR, MPI_COMM_WORLD);
+
+ /* For the centres we use a sum as MPI does not like bit-wise operations
+ on floating point numbers */
+ MPI_Allreduce(MPI_IN_PLACE, centres, 3 * nr_cells, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD);
+#endif
+
+ /* When writing a single file, only rank 0 writes the meta-data */
+ if ((distributed) || (!distributed && nodeID == 0)) {
+
+ /* Unit conversion if necessary */
+ const double factor = units_conversion_factor(
+ internal_units, snapshot_units, UNIT_CONV_LENGTH);
+ if (factor != 1.) {
+
+ /* Convert the cell centres */
+ for (int i = 0; i < nr_cells; ++i) {
+ centres[i * 3 + 0] *= factor;
+ centres[i * 3 + 1] *= factor;
+ centres[i * 3 + 2] *= factor;
+ }
+
+ /* Convert the cell widths */
+ cell_width[0] *= factor;
+ cell_width[1] *= factor;
+ cell_width[2] *= factor;
+ }
+
+ /* Write some meta-information first */
+ hid_t h_subgrp =
+ H5Gcreate(h_grp, "Meta-data", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_subgrp < 0) error("Error while creating meta-data sub-group");
+ io_write_attribute(h_subgrp, "nr_cells", INT, &nr_cells, 1);
+ io_write_attribute(h_subgrp, "size", DOUBLE, cell_width, 3);
+ io_write_attribute(h_subgrp, "dimension", INT, cdim, 3);
+ H5Gclose(h_subgrp);
+
+ /* Write the centres to the group */
+ hsize_t shape[2] = {(hsize_t)nr_cells, 3};
+ hid_t h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0) error("Error while creating data space for cell centres");
+ hid_t h_err = H5Sset_extent_simple(h_space, 2, shape, shape);
+ if (h_err < 0)
+ error("Error while changing shape of gas offsets data space.");
+ hid_t h_data = H5Dcreate(h_grp, "Centres", io_hdf5_type(DOUBLE), h_space,
+ H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace for gas offsets.");
+ h_err = H5Dwrite(h_data, io_hdf5_type(DOUBLE), h_space, H5S_ALL,
+ H5P_DEFAULT, centres);
+ if (h_err < 0) error("Error while writing centres.");
+ H5Dclose(h_data);
+ H5Sclose(h_space);
+
+ /* Group containing the offsets and counts for each particle type */
+ hid_t h_grp_offsets = H5Gcreate(h_grp, "OffsetsInFile", H5P_DEFAULT,
+ H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp_offsets < 0) error("Error while creating offsets sub-group");
+ hid_t h_grp_files =
+ H5Gcreate(h_grp, "Files", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp_files < 0) error("Error while creating filess sub-group");
+ hid_t h_grp_counts =
+ H5Gcreate(h_grp, "Counts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp_counts < 0) error("Error while creating counts sub-group");
+
+ if (global_counts[swift_type_gas] > 0 && num_fields[swift_type_gas] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType0", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_part, LONGLONG,
+ "PartType0", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_part, LONGLONG, "PartType0",
+ "counts");
+ }
+
+ if (global_counts[swift_type_dark_matter] > 0 &&
+ num_fields[swift_type_dark_matter] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType1", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_gpart, LONGLONG,
+ "PartType1", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_gpart, LONGLONG, "PartType1",
+ "counts");
+ }
+
+ if (global_counts[swift_type_dark_matter_background] > 0 &&
+ num_fields[swift_type_dark_matter_background] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType2", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_background_gpart, LONGLONG,
+ "PartType2", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_background_gpart, LONGLONG,
+ "PartType2", "counts");
+ }
+
+ if (global_counts[swift_type_sink] > 0 && num_fields[swift_type_sink] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType3", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_sink, LONGLONG,
+ "PartType3", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_sink, LONGLONG, "PartType3",
+ "counts");
+ }
+
+ if (global_counts[swift_type_stars] > 0 &&
+ num_fields[swift_type_stars] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType4", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_spart, LONGLONG,
+ "PartType4", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_spart, LONGLONG, "PartType4",
+ "counts");
+ }
+
+ if (global_counts[swift_type_black_hole] > 0 &&
+ num_fields[swift_type_black_hole] > 0) {
+ io_write_array(h_grp_files, nr_cells, files, INT, "PartType5", "files");
+ io_write_array(h_grp_offsets, nr_cells, offset_bpart, LONGLONG,
+ "PartType5", "offsets");
+ io_write_array(h_grp_counts, nr_cells, count_bpart, LONGLONG, "PartType5",
+ "counts");
+ }
+
+ H5Gclose(h_grp_offsets);
+ H5Gclose(h_grp_files);
+ H5Gclose(h_grp_counts);
+ }
+
+ /* Free everything we allocated */
+ free(centres);
+ free(files);
+ free(count_part);
+ free(count_gpart);
+ free(count_background_gpart);
+ free(count_spart);
+ free(count_bpart);
+ free(count_sink);
+ free(offset_part);
+ free(offset_gpart);
+ free(offset_background_gpart);
+ free(offset_spart);
+ free(offset_bpart);
+ free(offset_sink);
+}
+
+#endif /* HAVE_HDF5 */
diff --git a/src/common_io_copy.c b/src/common_io_copy.c
new file mode 100644
index 0000000000000000000000000000000000000000..e7e9f157726591a8e19c15c51df947e5617132ca
--- /dev/null
+++ b/src/common_io_copy.c
@@ -0,0 +1,753 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "common_io.h"
+
+/* Local includes. */
+#include "engine.h"
+#include "io_properties.h"
+#include "threadpool.h"
+
+/**
+ * @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 = (char*)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 xpart* restrict xparts = props.xparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ float* restrict temp_f = (float*)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, xparts + delta + i,
+ &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #part into a buffer of ints using a
+ * conversion function.
+ */
+void io_convert_part_i_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 xpart* restrict xparts = props.xparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ int* restrict temp_i = (int*)temp;
+ const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_part_i(e, parts + delta + i, xparts + delta + i,
+ &temp_i[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 xpart* restrict xparts = props.xparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ double* restrict temp_d = (double*)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, xparts + delta + i,
+ &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #part into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_part_l_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 xpart* restrict xparts = props.xparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_part_l(e, parts + delta + i, xparts + delta + i,
+ &temp_l[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 = (float*)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 ints using a
+ * conversion function.
+ */
+void io_convert_gpart_i_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? */
+ int* restrict temp_i = (int*)temp;
+ const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_gpart_i(e, gparts + delta + i, &temp_i[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 = (double*)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 Mapper function to copy #gpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_gpart_l_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? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_gpart_l(e, gparts + delta + i, &temp_l[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #spart into a buffer of floats using a
+ * conversion function.
+ */
+void io_convert_spart_f_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct spart* restrict sparts = props.sparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ float* restrict temp_f = (float*)temp;
+ const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_spart_f(e, sparts + delta + i, &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #spart into a buffer of ints using a
+ * conversion function.
+ */
+void io_convert_spart_i_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct spart* restrict sparts = props.sparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ int* restrict temp_i = (int*)temp;
+ const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_spart_i(e, sparts + delta + i, &temp_i[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #spart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_spart_d_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct spart* restrict sparts = props.sparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ double* restrict temp_d = (double*)temp;
+ const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_spart_d(e, sparts + delta + i, &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #spart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_spart_l_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct spart* restrict sparts = props.sparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_spart_l(e, sparts + delta + i, &temp_l[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of floats using a
+ * conversion function.
+ */
+void io_convert_bpart_f_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ float* restrict temp_f = (float*)temp;
+ const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_f(e, bparts + delta + i, &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of ints using a
+ * conversion function.
+ */
+void io_convert_bpart_i_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ int* restrict temp_i = (int*)temp;
+ const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_i(e, bparts + delta + i, &temp_i[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_bpart_d_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ double* restrict temp_d = (double*)temp;
+ const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_d(e, bparts + delta + i, &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #bpart into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_bpart_l_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct bpart* restrict bparts = props.bparts;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_bpart_l(e, bparts + delta + i, &temp_l[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of floats using a
+ * conversion function.
+ */
+void io_convert_sink_f_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct sink* restrict sinks = props.sinks;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ float* restrict temp_f = (float*)temp;
+ const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_sink_f(e, sinks + delta + i, &temp_f[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of ints using a
+ * conversion function.
+ */
+void io_convert_sink_i_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct sink* restrict sinks = props.sinks;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ int* restrict temp_i = (int*)temp;
+ const ptrdiff_t delta = (temp_i - props.start_temp_i) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_sink_i(e, sinks + delta + i, &temp_i[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_sink_d_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct sink* restrict sinks = props.sinks;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ double* restrict temp_d = (double*)temp;
+ const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_sink_d(e, sinks + delta + i, &temp_d[i * dim]);
+}
+
+/**
+ * @brief Mapper function to copy #sink into a buffer of doubles using a
+ * conversion function.
+ */
+void io_convert_sink_l_mapper(void* restrict temp, int N,
+ void* restrict extra_data) {
+
+ const struct io_props props = *((const struct io_props*)extra_data);
+ const struct sink* restrict sinks = props.sinks;
+ const struct engine* e = props.e;
+ const size_t dim = props.dimension;
+
+ /* How far are we with this chunk? */
+ long long* restrict temp_l = (long long*)temp;
+ const ptrdiff_t delta = (temp_l - props.start_temp_l) / dim;
+
+ for (int i = 0; i < N; i++)
+ props.convert_sink_l(e, sinks + delta + i, &temp_l[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 = (char*)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, threadpool_auto_chunk_size, (void*)&props);
+
+ } else { /* Converting particle to data */
+
+ if (props.convert_part_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)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,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_part_i != NULL) {
+
+ /* Prepare some parameters */
+ int* temp_i = (int*)temp;
+ props.start_temp_i = (int*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_part_i_mapper, temp_i, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_part_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)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,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_part_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_part_l_mapper, temp_l, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_gpart_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)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,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_gpart_i != NULL) {
+
+ /* Prepare some parameters */
+ int* temp_i = (int*)temp;
+ props.start_temp_i = (int*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_gpart_i_mapper, temp_i, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_gpart_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)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,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_gpart_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_gpart_l_mapper, temp_l, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_spart_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_spart_f_mapper, temp_f, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_spart_i != NULL) {
+
+ /* Prepare some parameters */
+ int* temp_i = (int*)temp;
+ props.start_temp_i = (int*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_spart_i_mapper, temp_i, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_spart_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_spart_d_mapper, temp_d, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_spart_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_spart_l_mapper, temp_l, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_sink_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_sink_f_mapper, temp_f, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_sink_i != NULL) {
+
+ /* Prepare some parameters */
+ int* temp_i = (int*)temp;
+ props.start_temp_i = (int*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_sink_i_mapper, temp_i, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_sink_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_sink_d_mapper, temp_d, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_sink_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_sink_l_mapper, temp_l, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_bpart_f != NULL) {
+
+ /* Prepare some parameters */
+ float* temp_f = (float*)temp;
+ props.start_temp_f = (float*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_f_mapper, temp_f, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_bpart_i != NULL) {
+
+ /* Prepare some parameters */
+ int* temp_i = (int*)temp;
+ props.start_temp_i = (int*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_i_mapper, temp_i, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_bpart_d != NULL) {
+
+ /* Prepare some parameters */
+ double* temp_d = (double*)temp;
+ props.start_temp_d = (double*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_d_mapper, temp_d, N, copySize,
+ threadpool_auto_chunk_size, (void*)&props);
+
+ } else if (props.convert_bpart_l != NULL) {
+
+ /* Prepare some parameters */
+ long long* temp_l = (long long*)temp;
+ props.start_temp_l = (long long*)temp;
+ props.e = e;
+
+ /* Copy the whole thing into a buffer */
+ threadpool_map((struct threadpool*)&e->threadpool,
+ io_convert_bpart_l_mapper, temp_l, N, copySize,
+ threadpool_auto_chunk_size, (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, (double*)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, (float*)temp,
+ IO_BUFFER_ALIGNMENT);
+ for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
+ }
+ }
+}
diff --git a/src/common_io_fields.c b/src/common_io_fields.c
new file mode 100644
index 0000000000000000000000000000000000000000..1eb9acfd75268ee83a7d865bd7c55821abc7f2bf
--- /dev/null
+++ b/src/common_io_fields.c
@@ -0,0 +1,389 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "common_io.h"
+
+/* Local includes. */
+#include "error.h"
+#include "units.h"
+
+/* I/O functions of each sub-module */
+#include "black_holes_io.h"
+#include "chemistry_io.h"
+#include "cooling_io.h"
+#include "fof_io.h"
+#include "gravity_io.h"
+#include "hydro_io.h"
+#include "rt_io.h"
+#include "sink_io.h"
+#include "star_formation_io.h"
+#include "stars_io.h"
+#include "tracers_io.h"
+#include "velociraptor_io.h"
+
+/* Some standard headers. */
+#include <string.h>
+
+/**
+ * @brief Return the particle type code of a select_output parameter
+ *
+ * @param name The name of the parameter under consideration.
+ *
+ * @return The (integer) particle type of the parameter.
+ */
+int io_get_param_ptype(const char* name) {
+
+ const int name_len = strlen(name);
+
+ for (int ptype = 0; ptype < swift_type_count; ptype++) {
+ const int ptype_name_len = strlen(part_type_names[ptype]);
+ if (name_len >= ptype_name_len &&
+ strcmp(&name[name_len - ptype_name_len], part_type_names[ptype]) == 0)
+ return ptype;
+ }
+
+ /* If we get here, we could not match the name, so something's gone wrong. */
+ error("Could not determine the particle type for parameter '%s'.", name);
+
+ /* We can never get here, but the compiler may complain if we don't return
+ * an int after promising to do so... */
+ return -1;
+}
+
+/**
+ * @brief Return the number and names of all output fields of a given ptype.
+ *
+ * @param ptype The index of the particle type under consideration.
+ * @param list An io_props list that will hold the individual fields.
+ * @param with_cosmology Use cosmological name variant?
+ * @param with_fof Include FoF related fields?
+ * @param with_stf Include STF related fields?
+ *
+ * @return The total number of fields that can be written for the ptype.
+ */
+int io_get_ptype_fields(const int ptype, struct io_props* list,
+ const int with_cosmology, const int with_fof,
+ const int with_stf) {
+
+ int num_fields = 0;
+
+ switch (ptype) {
+
+ case swift_type_gas:
+ hydro_write_particles(NULL, NULL, list, &num_fields);
+ num_fields += chemistry_write_particles(NULL, NULL, list + num_fields,
+ with_cosmology);
+ num_fields +=
+ cooling_write_particles(NULL, NULL, list + num_fields, NULL);
+ num_fields += tracers_write_particles(NULL, NULL, list + num_fields,
+ with_cosmology);
+ num_fields +=
+ star_formation_write_particles(NULL, NULL, list + num_fields);
+ if (with_fof)
+ num_fields += fof_write_parts(NULL, NULL, list + num_fields);
+ if (with_stf)
+ num_fields += velociraptor_write_parts(NULL, NULL, list + num_fields);
+ num_fields += rt_write_particles(NULL, list + num_fields);
+ break;
+
+ case swift_type_dark_matter:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ if (with_fof) num_fields += fof_write_gparts(NULL, list + num_fields);
+ if (with_stf)
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
+ case swift_type_dark_matter_background:
+ darkmatter_write_particles(NULL, list, &num_fields);
+ if (with_fof) num_fields += fof_write_gparts(NULL, list + num_fields);
+ if (with_stf)
+ num_fields += velociraptor_write_gparts(NULL, list + num_fields);
+ break;
+
+ case swift_type_stars:
+ stars_write_particles(NULL, list, &num_fields, with_cosmology);
+ num_fields += chemistry_write_sparticles(NULL, list + num_fields);
+ num_fields +=
+ tracers_write_sparticles(NULL, list + num_fields, with_cosmology);
+ num_fields += star_formation_write_sparticles(NULL, list + num_fields);
+ if (with_fof) num_fields += fof_write_sparts(NULL, list + num_fields);
+ if (with_stf)
+ num_fields += velociraptor_write_sparts(NULL, list + num_fields);
+ num_fields += rt_write_stars(NULL, list + num_fields);
+ break;
+
+ case swift_type_sink:
+ sink_write_particles(NULL, list, &num_fields, with_cosmology);
+ break;
+
+ case swift_type_black_hole:
+ black_holes_write_particles(NULL, list, &num_fields, with_cosmology);
+ num_fields += chemistry_write_bparticles(NULL, list + num_fields);
+ if (with_fof) num_fields += fof_write_bparts(NULL, list + num_fields);
+ if (with_stf)
+ num_fields += velociraptor_write_bparts(NULL, list + num_fields);
+ break;
+
+ default:
+ error("Particle Type %d not yet supported. Aborting", ptype);
+ }
+
+ return num_fields;
+}
+
+/**
+ * @brief Prepare the output option fields according to the user's choices and
+ * verify that they are valid.
+ *
+ * @param output_options The #output_options for this run
+ * @param with_cosmology Ran with cosmology?
+ * @param with_fof Are we running with on-the-fly Fof?
+ * @param with_stf Are we running with on-the-fly structure finder?
+ */
+void io_prepare_output_fields(struct output_options* output_options,
+ const int with_cosmology, const int with_fof,
+ const int with_stf) {
+
+ const int MAX_NUM_PTYPE_FIELDS = 100;
+
+ /* Parameter struct for the output options */
+ struct swift_params* params = output_options->select_output;
+
+ /* Get all possible outputs per particle type */
+ int ptype_num_fields_total[swift_type_count] = {0};
+ struct io_props field_list[swift_type_count][MAX_NUM_PTYPE_FIELDS];
+
+ for (int ptype = 0; ptype < swift_type_count; ptype++)
+ ptype_num_fields_total[ptype] = io_get_ptype_fields(
+ ptype, field_list[ptype], with_cosmology, with_fof, with_stf);
+
+ /* Check for whether we have a `Default` section */
+ int have_default = 0;
+
+ /* Loop over each section, i.e. different class of output */
+ for (int section_id = 0; section_id < params->sectionCount; section_id++) {
+
+ /* Get the name of current (selection) section, without a trailing colon */
+ char section_name[FIELD_BUFFER_SIZE];
+ strcpy(section_name, params->section[section_id].name);
+ section_name[strlen(section_name) - 1] = 0;
+
+ /* Is this the `Default` section? */
+ if (strcmp(section_name, select_output_header_default_name) == 0)
+ have_default = 1;
+
+ /* How many fields should each ptype write by default? */
+ int ptype_num_fields_to_write[swift_type_count];
+
+ /* What is the default writing status for each ptype (on/off)? */
+ int ptype_default_write_status[swift_type_count];
+
+ /* Initialise section-specific writing counters for each particle type.
+ * If default is 'write', then we start from the total to deduct any fields
+ * that are switched off. If the default is 'off', we have to start from
+ * zero and then count upwards for each field that is switched back on. */
+ for (int ptype = 0; ptype < swift_type_count; ptype++) {
+
+ /* Internally also verifies that the default level is allowed */
+ const enum lossy_compression_schemes compression_level_current_default =
+ output_options_get_ptype_default_compression(params, section_name,
+ (enum part_type)ptype);
+
+ if (compression_level_current_default == compression_do_not_write) {
+ ptype_default_write_status[ptype] = 0;
+ ptype_num_fields_to_write[ptype] = 0;
+ } else {
+ ptype_default_write_status[ptype] = 1;
+ ptype_num_fields_to_write[ptype] = ptype_num_fields_total[ptype];
+ }
+
+ } /* ends loop over particle types */
+
+ /* Loop over each parameter */
+ for (int param_id = 0; param_id < params->paramCount; param_id++) {
+
+ /* Full name of the parameter to check */
+ const char* param_name = params->data[param_id].name;
+
+ /* Check whether the file still contains the old, now inappropriate
+ * 'SelectOutput' section */
+ if (strstr(param_name, "SelectOutput:") != NULL) {
+ error(
+ "Output selection files no longer require the use of top level "
+ "SelectOutput; see the documentation for changes.");
+ }
+
+ /* Skip if the parameter belongs to another output class or is a
+ * 'Standard' parameter */
+ if (strstr(param_name, section_name) == NULL) continue;
+ if (strstr(param_name, ":Standard_") != NULL) continue;
+
+ /* Get the particle type for current parameter
+ * (raises an error if it could not determine it) */
+ const int param_ptype = io_get_param_ptype(param_name);
+
+ /* Issue a warning if this parameter does not pertain to any of the
+ * known fields from this ptype. */
+ int field_id = 0;
+ char field_name[PARSER_MAX_LINE_SIZE];
+ for (field_id = 0; field_id < ptype_num_fields_total[param_ptype];
+ field_id++) {
+
+ sprintf(field_name, "%s:%.*s_%s", section_name, FIELD_BUFFER_SIZE,
+ field_list[param_ptype][field_id].name,
+ part_type_names[param_ptype]);
+
+ if (strcmp(param_name, field_name) == 0) break;
+ }
+
+ int param_is_known = 0; /* Update below if it is a known one */
+ if (field_id < ptype_num_fields_total[param_ptype])
+ param_is_known = 1;
+ else
+ message(
+ "WARNING: Trying to change behaviour of field '%s' (read from "
+ "'%s') that does not exist. This may be because you are not "
+ "running with all of the physics that you compiled the code with.",
+ param_name, params->fileName);
+
+ /* Perform a correctness check on the _value_ of the parameter */
+ char param_value[FIELD_BUFFER_SIZE];
+ parser_get_param_string(params, param_name, param_value);
+
+ int value_id = 0;
+ for (value_id = 0; value_id < compression_level_count; value_id++)
+ if (strcmp(param_value, lossy_compression_schemes_names[value_id]) == 0)
+ break;
+
+ if (value_id == compression_level_count)
+ error("Choice of output selection parameter %s ('%s') is invalid.",
+ param_name, param_value);
+
+ /* Adjust number of fields to be written for param_ptype, if this field's
+ * status is different from default and it is a known one. */
+ if (param_is_known) {
+ const int is_on =
+ strcmp(param_value,
+ lossy_compression_schemes_names[compression_do_not_write]) !=
+ 0;
+
+ if (is_on && !ptype_default_write_status[param_ptype]) {
+ /* Particle should be written even though default is off:
+ * increase field count */
+ ptype_num_fields_to_write[param_ptype] += 1;
+ }
+ if (!is_on && ptype_default_write_status[param_ptype]) {
+ /* Particle should not be written, even though default is on:
+ * decrease field count */
+ ptype_num_fields_to_write[param_ptype] -= 1;
+ }
+ }
+ } /* ends loop over parameters */
+
+ /* Second loop over ptypes, to write out total number of fields to write */
+ for (int ptype = 0; ptype < swift_type_count; ptype++) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Sanity check: is the number of fields to write non-negative? */
+ if (ptype_num_fields_to_write[ptype] < 0)
+ error(
+ "We seem to have subtracted too many fields for particle "
+ "type %d in output class %s (total to write is %d)",
+ ptype, section_name, ptype_num_fields_to_write[ptype]);
+#endif
+ output_options->num_fields_to_write[section_id][ptype] =
+ ptype_num_fields_to_write[ptype];
+ }
+ } /* Ends loop over sections, for different output classes */
+
+ /* Add field numbers for (possible) implicit `Default` output class */
+ if (!have_default) {
+ const int default_id = output_options->select_output->sectionCount;
+ for (int ptype = 0; ptype < swift_type_count; ptype++)
+ output_options->num_fields_to_write[default_id][ptype] =
+ ptype_num_fields_total[ptype];
+ }
+}
+
+/**
+ * @brief Write the output field parameters file
+ *
+ * @param filename The file to write.
+ * @param with_cosmology Use cosmological name variant?
+ */
+void io_write_output_field_parameter(const char* filename, int with_cosmology) {
+
+ FILE* file = fopen(filename, "w");
+ if (file == NULL) error("Error opening file '%s'", filename);
+
+ /* Create a fake unit system for the snapshots */
+ struct unit_system snapshot_units;
+ units_init_cgs(&snapshot_units);
+
+ /* Loop over all particle types */
+ fprintf(file, "Default:\n");
+ for (int ptype = 0; ptype < swift_type_count; ptype++) {
+
+ struct io_props list[100];
+ int num_fields = io_get_ptype_fields(ptype, list, with_cosmology,
+ /*with_fof=*/1, /*with_stf=*/1);
+
+ if (num_fields == 0) continue;
+
+ /* Output a header for that particle type */
+ fprintf(file, " # Particle Type %s\n", part_type_names[ptype]);
+
+ /* Write all the fields of this particle type */
+ for (int i = 0; i < num_fields; ++i) {
+
+ char unit_buffer[FIELD_BUFFER_SIZE] = {0};
+ units_cgs_conversion_string(unit_buffer, &snapshot_units, list[i].units,
+ list[i].scale_factor_exponent);
+
+ /* Need to buffer with a maximal size - otherwise we can't read in again
+ * because comments are too long */
+ char comment_write_buffer[PARSER_MAX_LINE_SIZE / 2];
+
+ sprintf(comment_write_buffer, "%.*s", PARSER_MAX_LINE_SIZE / 2 - 1,
+ list[i].description);
+
+ /* If our string is too long, replace the last few characters (before
+ * \0) with ... for 'fancy printing' */
+ if (strlen(comment_write_buffer) > PARSER_MAX_LINE_SIZE / 2 - 3) {
+ strcpy(&comment_write_buffer[PARSER_MAX_LINE_SIZE / 2 - 4], "...");
+ }
+
+ fprintf(file, " %s_%s: %s # %s : %s\n", list[i].name,
+ part_type_names[ptype], "on", comment_write_buffer, unit_buffer);
+ }
+
+ fprintf(file, "\n");
+ }
+
+ fclose(file);
+
+ printf(
+ "List of valid ouput fields for the particle in snapshots dumped in "
+ "'%s'.\n",
+ filename);
+}
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 4e98c24397d51da2e1a43899aa40ab22f43db96d..55b345912da6b57060bf587179089afcf3c88f9a 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -1046,7 +1046,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
}
/* If we are remapping ParticleIDs later, start by setting them to 1. */
- if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+ if (remap_ids) io_set_ids_to_one(*gparts, *Ngparts);
if (!dry_run && with_gravity) {
diff --git a/src/serial_io.c b/src/serial_io.c
index b734d5fa7082b059537887f3260ceafa7c1a22a6..eb801c8f10a21f1336e6dc13897f57925a0905bb 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -851,7 +851,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
}
/* If we are remapping ParticleIDs later, start by setting them to 1. */
- if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+ if (remap_ids) io_set_ids_to_one(*gparts, *Ngparts);
/* Duplicate the parts for gravity */
if (!dry_run && with_gravity) {
diff --git a/src/single_io.c b/src/single_io.c
index 3765257f3914a0610763ad611ac2d8aef99c88cb..e0912ba45e53945a3ee98d9d945343df77292b53 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -708,7 +708,7 @@ void read_ic_single(const char* fileName,
}
/* If we are remapping ParticleIDs later, start by setting them to 1. */
- if (remap_ids) set_ids_to_one(*gparts, *Ngparts);
+ if (remap_ids) io_set_ids_to_one(*gparts, *Ngparts);
/* Duplicate the parts for gravity */
if (!dry_run && with_gravity) {