diff --git a/src/cell.c b/src/cell.c
index 16acc4c4602a088b8b33cd558c0b941ee9429c23..fceaed8ac83a7754114179d35c7c8f91ad6cf262 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -4524,9 +4524,7 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
/* Check that we are actually going to move forward. */
- if (ti_current < ti_old_part)
- error("Attempt to drift to the past ti_current=%lld < ti_old_part=%lld",
- ti_current, ti_old_part);
+ if (ti_current < ti_old_part) error("Attempt to drift to the past");
#endif
/* Early abort? */
diff --git a/src/engine.c b/src/engine.c
index 758cd374efba1db8b0ccccbffcd8f254d3485b51..48aa501d165eaff1514ab2dcdfc2f346b35cc131 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2902,8 +2902,11 @@ void engine_check_for_dumps(struct engine *e) {
break;
case output_los:
+
+ /* Compute the LoS */
do_line_of_sight(e);
+ /* Move on */
engine_compute_next_los_time(e);
break;
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index b17f0a25a188708d7e299f086124242fb8f45f4f..f5d2557bc2a942e6d8413883f3f8bc09206e29a0 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -43,6 +43,61 @@
#include "tracers_io.h"
#include "velociraptor_io.h"
+/**
+ * @brief Will the line of sight intersect a given cell?
+ *
+ * Also return 0 if the cell is empty.
+ *
+ * @param c The top level cell.
+ * @param los The line of sight structure.
+ */
+static INLINE int does_los_intersect(const struct cell *c,
+ const struct line_of_sight *los) {
+
+ /* Empty cell? */
+ if (c->hydro.count == 0) return 0;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->hydro.h_max <= 0.) error("Invalid h_max for does_los_intersect");
+#endif
+
+ /* Distance from LOS to left and bottom cell edge. */
+ const double cx_min = c->loc[los->xaxis];
+ const double cy_min = c->loc[los->yaxis];
+
+ /* Distance from LOS to right and top cell edge. */
+ const double cx_max = c->loc[los->xaxis] + c->width[los->xaxis];
+ const double cy_max = c->loc[los->yaxis] + c->width[los->yaxis];
+
+ /* Maximum smoothing length of a part in this top level cell. */
+ const double hsml = c->hydro.h_max * kernel_gamma;
+ const double hsml2 = hsml * hsml;
+
+ double dx, dy;
+
+ if (los->periodic) {
+ dx = min(fabs(nearest(cx_min - los->Xpos, los->dim[los->xaxis])),
+ fabs(nearest(cx_max - los->Xpos, los->dim[los->xaxis])));
+ dy = min(fabs(nearest(cy_min - los->Ypos, los->dim[los->yaxis])),
+ fabs(nearest(cy_max - los->Ypos, los->dim[los->yaxis])));
+ } else {
+ dx = fabs(cx_max - los->Xpos);
+ dy = fabs(cy_max - los->Ypos);
+ }
+
+ /* Is sightline directly within this top level cell? */
+ if (dx < (1.01 * c->width[los->xaxis]) / 2. &&
+ dy < (1.01 * c->width[los->yaxis]) / 2.) {
+ return 1;
+ /* Could a part from this top level cell smooth into the sightline? */
+ } else if (dx * dx + dy * dy < hsml2) {
+ return 1;
+ /* Don't need to work with this top level cell. */
+ } else {
+ return 0;
+ }
+}
+
/**
* @brief Reads the LOS properties from the param file.
*
@@ -50,8 +105,9 @@
* @param los_params Sightline parameters to save into.
* @param params Swift params to read from.
*/
-void los_init(double dim[3], struct los_props *los_params,
+void los_init(const double dim[3], struct los_props *los_params,
struct swift_params *params) {
+
/* How many line of sights in each plane. */
los_params->num_along_z =
parser_get_opt_param_int(params, "LineOfSight:num_along_z", 0);
@@ -101,6 +157,30 @@ void los_init(double dim[3], struct los_props *los_params,
los_params->range_when_shooting_down_axis[2]);
}
+/**
+ * @brief Create a #line_of_sight object from its attributes
+ */
+void create_sightline(const double Xpos, const double Ypos,
+ enum los_direction xaxis, enum los_direction yaxis,
+ enum los_direction zaxis, const int periodic,
+ const double dim[3], struct line_of_sight *los,
+ const double range_when_shooting_down_axis[2]) {
+ los->Xpos = Xpos;
+ los->Ypos = Ypos;
+ los->particles_in_los_local = 0;
+ los->particles_in_los_total = 0;
+ los->xaxis = xaxis;
+ los->yaxis = yaxis;
+ los->zaxis = zaxis;
+ los->periodic = periodic;
+ los->dim[0] = dim[0];
+ los->dim[1] = dim[1];
+ los->dim[2] = dim[2];
+ los->num_intersecting_top_level_cells = 0;
+ los->range_when_shooting_down_axis[0] = range_when_shooting_down_axis[0];
+ los->range_when_shooting_down_axis[1] = range_when_shooting_down_axis[1];
+}
+
/**
* @brief Generates random sightline positions.
*
@@ -116,48 +196,55 @@ void generate_sightlines(struct line_of_sight *Los,
/* Keep track of number of sightlines. */
int count = 0;
- double Xpos, Ypos;
-
/* Sightlines in XY plane, shoots down Z. */
for (int i = 0; i < params->num_along_z; i++) {
- Xpos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
- params->allowed_losrange_x[0];
- Ypos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
- params->allowed_losrange_y[0];
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
+ params->allowed_losrange_x[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
+ params->allowed_losrange_y[0];
+
create_sightline(Xpos, Ypos, simulation_x_axis, simulation_y_axis,
simulation_z_axis, periodic, dim, &Los[count],
params->range_when_shooting_down_axis[simulation_z_axis]);
- count += 1;
+ count++;
}
/* Sightlines in YZ plane, shoots down X. */
for (int i = 0; i < params->num_along_x; i++) {
- Xpos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
- params->allowed_losrange_y[0];
- Ypos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
- params->allowed_losrange_z[0];
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_y[1] - params->allowed_losrange_y[0])) +
+ params->allowed_losrange_y[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
+ params->allowed_losrange_z[0];
+
create_sightline(Xpos, Ypos, simulation_y_axis, simulation_z_axis,
simulation_x_axis, periodic, dim, &Los[count],
params->range_when_shooting_down_axis[simulation_x_axis]);
- count += 1;
+ count++;
}
/* Sightlines in XZ plane, shoots down Y. */
for (int i = 0; i < params->num_along_y; i++) {
- Xpos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
- params->allowed_losrange_x[0];
- Ypos = ((float)rand() / (float)(RAND_MAX) *
- (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
- params->allowed_losrange_z[0];
+ double Xpos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_x[1] - params->allowed_losrange_x[0])) +
+ params->allowed_losrange_x[0];
+ double Ypos =
+ ((float)rand() / (float)(RAND_MAX) *
+ (params->allowed_losrange_z[1] - params->allowed_losrange_z[0])) +
+ params->allowed_losrange_z[0];
+
create_sightline(Xpos, Ypos, simulation_x_axis, simulation_z_axis,
simulation_y_axis, periodic, dim, &Los[count],
params->range_when_shooting_down_axis[simulation_y_axis]);
- count += 1;
+ count++;
}
/* Make sure we made the correct ammount */
@@ -165,27 +252,6 @@ void generate_sightlines(struct line_of_sight *Los,
error("Could not make the right number of sightlines");
}
-void create_sightline(const double Xpos, const double Ypos,
- enum los_direction xaxis, enum los_direction yaxis,
- enum los_direction zaxis, const int periodic,
- const double dim[3], struct line_of_sight *los,
- const double range_when_shooting_down_axis[2]) {
- los->Xpos = Xpos;
- los->Ypos = Ypos;
- los->particles_in_los_local = 0;
- los->particles_in_los_total = 0;
- los->xaxis = xaxis;
- los->yaxis = yaxis;
- los->zaxis = zaxis;
- los->periodic = periodic;
- los->dim[0] = dim[0];
- los->dim[1] = dim[1];
- los->dim[2] = dim[2];
- los->num_intersecting_top_level_cells = 0;
- los->range_when_shooting_down_axis[0] = range_when_shooting_down_axis[0];
- los->range_when_shooting_down_axis[1] = range_when_shooting_down_axis[1];
-}
-
/**
* @brief Print line_of_sight information.
*
@@ -194,88 +260,372 @@ void create_sightline(const double Xpos, const double Ypos,
void print_los_info(const struct line_of_sight *Los, const int i) {
message(
- "[LOS %i] Xpos:%g Ypos:%g parts_in_los:%li "
+ "[LOS %i] Xpos:%g Ypos:%g parts_in_los:%i "
"num_intersecting_top_level_cells:%i",
i, Los[i].Xpos, Los[i].Ypos, Los[i].particles_in_los_total,
Los[i].num_intersecting_top_level_cells);
}
/**
- * @brief Loop over each part to see which ones intersect the LOS.
+ * @brief Writes dataset for a given part attribute.
*
- * @param map_data The parts.
- * @param count The number of parts.
- * @param extra_data The line_of_sight structure for this LOS.
+ * @param p io_props dataset for this attribute.
+ * @param N number of parts in this line of sight.
+ * @param j Line of sight ID.
+ * @param e The engine.
+ * @param grp HDF5 group to write to.
*/
-void los_first_loop_mapper(void *restrict map_data, int count,
- void *restrict extra_data) {
+void write_los_hdf5_dataset(const struct io_props props, const size_t N,
+ const int j, const struct engine *e, hid_t grp) {
- size_t los_particle_count = 0;
- double dx, dy, dx2, dy2, hsml, hsml2;
- struct line_of_sight *restrict LOS_list = (struct line_of_sight *)extra_data;
- struct part *restrict parts = (struct part *)map_data;
+ /* Create data space */
+ const hid_t h_space = H5Screate(H5S_SIMPLE);
+ if (h_space < 0)
+ error("Error while creating data space for field '%s'.", props.name);
- /* Loop over each part to find those in LOS. */
- for (int i = 0; i < count; i++) {
+ int rank = 0;
+ hsize_t shape[2];
+ hsize_t chunk_shape[2];
+ if (props.dimension > 1) {
+ rank = 2;
+ shape[0] = N;
+ shape[1] = props.dimension;
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
+ chunk_shape[1] = props.dimension;
+ } else {
+ rank = 1;
+ shape[0] = N;
+ shape[1] = 0;
+ chunk_shape[0] = 1 << 20; /* Just a guess...*/
+ chunk_shape[1] = 0;
+ }
- /* Don't consider inhibited parts. */
- if (parts[i].time_bin == time_bin_inhibited) continue;
+ /* Make sure the chunks are not larger than the dataset */
+ if (chunk_shape[0] > N) chunk_shape[0] = N;
- /* Don't consider part if outwith allowed z-range. */
- if (parts[i].x[LOS_list->zaxis] <
- LOS_list->range_when_shooting_down_axis[0] ||
- parts[i].x[LOS_list->zaxis] >
- LOS_list->range_when_shooting_down_axis[1])
- continue;
+ /* Change shape of data space */
+ hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
+ if (h_err < 0)
+ error("Error while changing data space shape for field '%s'.", props.name);
- /* Distance from this part to LOS along x dim. */
- dx = parts[i].x[LOS_list->xaxis] - LOS_list->Xpos;
+ /* Dataset properties */
+ const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
- /* Periodic wrap. */
- if (LOS_list->periodic) dx = nearest(dx, LOS_list->dim[LOS_list->xaxis]);
+ /* Set chunk size */
+ h_err = H5Pset_chunk(h_prop, rank, chunk_shape);
+ if (h_err < 0)
+ error("Error while setting chunk size (%llu, %llu) for field '%s'.",
+ chunk_shape[0], chunk_shape[1], props.name);
- /* Square. */
- dx2 = dx * dx;
+ /* Impose check-sum to verify data corruption */
+ h_err = H5Pset_fletcher32(h_prop);
+ if (h_err < 0)
+ error("Error while setting checksum options for field '%s'.", props.name);
- /* Smoothing length of this part. */
- hsml = parts[i].h * kernel_gamma;
- hsml2 = hsml * hsml;
+ /* Impose data compression */
+ if (e->snapshot_compression > 0) {
+ h_err = H5Pset_shuffle(h_prop);
+ if (h_err < 0)
+ error("Error while setting shuffling options for field '%s'.",
+ props.name);
- /* Does this particle fall into our LOS? */
- if (dx2 < hsml2) {
+ h_err = H5Pset_deflate(h_prop, e->snapshot_compression);
+ if (h_err < 0)
+ error("Error while setting compression options for field '%s'.",
+ props.name);
+ }
- /* Distance from this part to LOS along y dim. */
- dy = parts[i].x[LOS_list->yaxis] - LOS_list->Ypos;
+ /* Allocate temporary buffer */
+ const size_t num_elements = N * props.dimension;
+ const size_t typeSize = io_sizeof_type(props.type);
+ void *temp = NULL;
+ if (swift_memalign("writebuff", (void **)&temp, IO_BUFFER_ALIGNMENT,
+ num_elements * typeSize) != 0)
+ error("Unable to allocate temporary i/o buffer");
- /* Periodic wrap. */
- if (LOS_list->periodic) dy = nearest(dy, LOS_list->dim[LOS_list->yaxis]);
+ /* Copy particle data to temp buffer */
+ io_copy_temp_buffer(temp, e, props, N, e->internal_units, e->snapshot_units);
- /* Square. */
- dy2 = dy * dy;
+ /* Create dataset */
+ char att_name[200];
+ sprintf(att_name, "/LOS_%04i/%s", j, props.name);
+ const hid_t h_data = H5Dcreate(grp, att_name, io_hdf5_type(props.type),
+ h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
+ if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
- /* Does this part still fall into our LOS? */
- if (dy2 < hsml2) {
+ /* Write dataset */
+ herr_t status = H5Dwrite(h_data, io_hdf5_type(props.type), H5S_ALL, H5S_ALL,
+ H5P_DEFAULT, temp);
+ if (status < 0) error("Error while writing data array '%s'.", props.name);
- /* 2D distance to LOS. */
- if (dx2 + dy2 <= hsml2) {
+ /* Write unit conversion factors for this data set */
+ char buffer[FIELD_BUFFER_SIZE] = {0};
+ units_cgs_conversion_string(buffer, e->snapshot_units, props.units,
+ props.scale_factor_exponent);
+ float baseUnitsExp[5];
+ units_get_base_unit_exponents_array(baseUnitsExp, props.units);
+ io_write_attribute_f(h_data, "U_M exponent", baseUnitsExp[UNIT_MASS]);
+ io_write_attribute_f(h_data, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
+ io_write_attribute_f(h_data, "U_t exponent", baseUnitsExp[UNIT_TIME]);
+ io_write_attribute_f(h_data, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
+ io_write_attribute_f(h_data, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
+ io_write_attribute_f(h_data, "h-scale exponent", 0.f);
+ io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
+ io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
- /* We've found one. */
- los_particle_count++;
- }
- }
- }
- } /* End of loop over all parts */
+ /* Write the actual number this conversion factor corresponds to */
+ const double factor =
+ units_cgs_conversion_factor(e->snapshot_units, props.units);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to CGS (not including cosmological corrections)",
+ factor);
+ io_write_attribute_d(
+ h_data,
+ "Conversion factor to physical CGS (including cosmological corrections)",
+ factor * pow(e->cosmology->a, props.scale_factor_exponent));
- atomic_add(&LOS_list->particles_in_los_local, los_particle_count);
+#ifdef SWIFT_DEBUG_CHECKS
+ if (strlen(props.description) == 0)
+ error("Invalid (empty) description of the field '%s'", props.name);
+#endif
+
+ /* Write the full description */
+ io_write_attribute_s(h_data, "Description", props.description);
+
+ /* Free and close everything */
+ swift_free("writebuff", temp);
+ H5Pclose(h_prop);
+ H5Dclose(h_data);
+ H5Sclose(h_space);
}
/**
- * @brief Find all top level cells that a LOS will intersect.
- *
- * This includes the top level cells the LOS directly passes through
- * and the neighbouring top level cells whose parts could smooth into the LOS.
+ * @brief Write parts in LOS to HDF5 file.
*
- * @param e The engine.
+ * @param grp HDF5 group of this LOS.
+ * @param j LOS ID.
+ * @param N number of parts in this line of sight.
+ * @param parts the list of parts in this LOS.
+ * @param e The engine.
+ * @param xparts the list of xparts in this LOS.
+ */
+void write_los_hdf5_datasets(hid_t grp, const int j, const size_t N,
+ const struct part *parts, const struct engine *e,
+ const struct xpart *xparts) {
+
+ /* What kind of run are we working with? */
+ struct swift_params *params = e->parameter_file;
+ const int with_cosmology = e->policy & engine_policy_cosmology;
+ const int with_cooling = e->policy & engine_policy_cooling;
+ const int with_temperature = e->policy & engine_policy_temperature;
+ const int with_fof = e->policy & engine_policy_fof;
+#ifdef HAVE_VELOCIRAPTOR
+ const int with_stf = (e->policy & engine_policy_structure_finding) &&
+ (e->s->gpart_group_data != NULL);
+#else
+ const int with_stf = 0;
+#endif
+
+ int num_fields = 0;
+ struct io_props list[100];
+
+ /* Find all the gas output fields */
+ hydro_write_particles(parts, xparts, list, &num_fields);
+ num_fields += chemistry_write_particles(parts, list + num_fields);
+ if (with_cooling || with_temperature) {
+ num_fields += cooling_write_particles(parts, xparts, list + num_fields,
+ e->cooling_func);
+ }
+ if (with_fof) {
+ num_fields += fof_write_parts(parts, xparts, list + num_fields);
+ }
+ if (with_stf) {
+ num_fields += velociraptor_write_parts(parts, xparts, list + num_fields);
+ }
+ num_fields +=
+ tracers_write_particles(parts, xparts, list + num_fields, with_cosmology);
+ num_fields +=
+ star_formation_write_particles(parts, xparts, list + num_fields);
+
+ /* Loop over each output field */
+ for (int i = 0; i < num_fields; i++) {
+
+ /* Did the user cancel this field? */
+ char field[PARSER_MAX_LINE_SIZE];
+ sprintf(field, "SelectOutputLOS:%.*s", FIELD_BUFFER_SIZE, list[i].name);
+ int should_write = parser_get_opt_param_int(params, field, 1);
+
+ /* Write (if selected) */
+ if (should_write) write_los_hdf5_dataset(list[i], N, j, e, grp);
+ }
+}
+
+/**
+ * @brief Writes HDF5 headers and information groups for this line of sight.
+ *
+ * @param h_file HDF5 file reference.
+ * @param e The engine.
+ * @param LOS_params The line of sight params.
+ */
+void write_hdf5_header(hid_t h_file, const struct engine *e,
+ const struct los_props *LOS_params,
+ const size_t total_num_parts_in_los) {
+
+ /* Open header to write simulation properties */
+ hid_t h_grp =
+ H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating file header\n");
+
+ /* Convert basic output information to snapshot units */
+ const double factor_time = units_conversion_factor(
+ e->internal_units, e->snapshot_units, UNIT_CONV_TIME);
+ const double factor_length = units_conversion_factor(
+ e->internal_units, e->snapshot_units, UNIT_CONV_LENGTH);
+ const double dblTime = e->time * factor_time;
+ const double dim[3] = {e->s->dim[0] * factor_length,
+ e->s->dim[1] * factor_length,
+ e->s->dim[2] * factor_length};
+
+ io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
+ io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
+ const int dimension = (int)hydro_dimension;
+ io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
+ io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
+ io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
+ io_write_attribute_s(h_grp, "Code", "SWIFT");
+ io_write_attribute_s(h_grp, "RunName", e->run_name);
+ const int num_files_per_snapshot = 1;
+ io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &num_files_per_snapshot,
+ 1);
+ io_write_attribute_s(h_grp, "OutputType", "LineOfSight");
+
+ /* GADGET-2 legacy values */
+ /* Number of particles of each type */
+ long long N_total[swift_type_count] = {0};
+ N_total[0] = total_num_parts_in_los;
+ unsigned int numParticles[swift_type_count] = {0};
+ unsigned int numParticlesHighWord[swift_type_count] = {0};
+ for (int ptype = 0; ptype < swift_type_count; ++ptype) {
+ numParticles[ptype] = (unsigned int)N_total[ptype];
+ numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
+ }
+ io_write_attribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total,
+ swift_type_count);
+ io_write_attribute(h_grp, "NumPart_Total", UINT, numParticles,
+ swift_type_count);
+ io_write_attribute(h_grp, "NumPart_Total_HighWord", UINT,
+ numParticlesHighWord, swift_type_count);
+
+ /* Store the time at which the snapshot was written */
+ time_t tm = time(NULL);
+ struct tm *timeinfo = localtime(&tm);
+ char snapshot_date[64];
+ strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
+ io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
+
+ /* Close group */
+ H5Gclose(h_grp);
+
+ io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units);
+
+ /* Print the LOS properties */
+ h_grp = H5Gcreate(h_file, "/LineOfSightParameters", H5P_DEFAULT, H5P_DEFAULT,
+ H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating LOS group");
+
+ /* Record this LOS attributes */
+ const int num_los_per_axis[3] = {LOS_params->num_along_x,
+ LOS_params->num_along_y,
+ LOS_params->num_along_z};
+ io_write_attribute(h_grp, "NumLineOfSight_PerAxis", INT, num_los_per_axis, 3);
+ io_write_attribute(h_grp, "NumLineOfSight_Total", INT, &LOS_params->num_tot,
+ 1);
+ io_write_attribute(h_grp, "AllowedLOSRangeX", DOUBLE,
+ LOS_params->allowed_losrange_x, 2);
+ io_write_attribute(h_grp, "AllowedLOSRangeY", DOUBLE,
+ LOS_params->allowed_losrange_y, 2);
+ io_write_attribute(h_grp, "AllowedLOSRangeZ", DOUBLE,
+ LOS_params->allowed_losrange_z, 2);
+ H5Gclose(h_grp);
+}
+
+/**
+ * @brief Loop over each part to see which ones intersect the LOS.
+ *
+ * @param map_data The parts.
+ * @param count The number of parts.
+ * @param extra_data The line_of_sight structure for this LOS.
+ */
+void los_first_loop_mapper(void *restrict map_data, int count,
+ void *restrict extra_data) {
+
+ struct line_of_sight *LOS_list = (struct line_of_sight *)extra_data;
+ const struct part *parts = (struct part *)map_data;
+
+ size_t los_particle_count = 0;
+
+ /* Loop over each part to find those in LOS. */
+ for (int i = 0; i < count; i++) {
+
+ /* Don't consider inhibited parts. */
+ if (parts[i].time_bin == time_bin_inhibited) continue;
+
+ /* Don't consider part if outwith allowed z-range. */
+ if (parts[i].x[LOS_list->zaxis] <
+ LOS_list->range_when_shooting_down_axis[0] ||
+ parts[i].x[LOS_list->zaxis] >
+ LOS_list->range_when_shooting_down_axis[1])
+ continue;
+
+ /* Distance from this part to LOS along x dim. */
+ double dx = parts[i].x[LOS_list->xaxis] - LOS_list->Xpos;
+
+ /* Periodic wrap. */
+ if (LOS_list->periodic) dx = nearest(dx, LOS_list->dim[LOS_list->xaxis]);
+
+ /* Square. */
+ const double dx2 = dx * dx;
+
+ /* Smoothing length of this part. */
+ const double hsml = parts[i].h * kernel_gamma;
+ const double hsml2 = hsml * hsml;
+
+ /* Does this particle fall into our LOS? */
+ if (dx2 < hsml2) {
+
+ /* Distance from this part to LOS along y dim. */
+ double dy = parts[i].x[LOS_list->yaxis] - LOS_list->Ypos;
+
+ /* Periodic wrap. */
+ if (LOS_list->periodic) dy = nearest(dy, LOS_list->dim[LOS_list->yaxis]);
+
+ /* Square. */
+ const double dy2 = dy * dy;
+
+ /* Does this part still fall into our LOS? */
+ if (dy2 < hsml2) {
+
+ /* 2D distance to LOS. */
+ if (dx2 + dy2 <= hsml2) {
+
+ /* We've found one. */
+ los_particle_count++;
+ }
+ }
+ }
+ } /* End of loop over all parts */
+
+ atomic_add(&LOS_list->particles_in_los_local, los_particle_count);
+}
+
+/**
+ * @brief Find all top level cells that a LOS will intersect.
+ *
+ * This includes the top level cells the LOS directly passes through
+ * and the neighbouring top level cells whose parts could smooth into the LOS.
+ *
+ * @param e The engine.
* @param los The line_of_sight structure.
*/
void find_intersecting_top_level_cells(
@@ -309,91 +659,42 @@ void find_intersecting_top_level_cells(
}
/**
- * @brief Will the line of sight intersect a given top level cell?
+ * @brief Main work function for computing line of sights.
*
- * @param c The top level cell.
- * @param los The line of sight structure.
+ * 1) Construct N random line of sight positions.
+ * 2) Loop over each line of sight.
+ * - 2.1) Find which top level cells sightline intersects.
+ * - 2.2) Loop over each part in these top level cells to see which intersect
+ * sightline.
+ * - 2.3) Use this count to construct a LOS parts/xparts array.
+ * - 2.4) Loop over each part and extract those in sightline to new array.
+ * - 2.5) Save sightline parts to HDF5 file.
+ *
+ * @param e The engine.
*/
-int does_los_intersect(const struct cell *c, const struct line_of_sight *los) {
+void do_line_of_sight(struct engine *e) {
- /* Empty cell? */
- if (c->hydro.count == 0) return 0;
+ /* Start counting. */
+ const ticks tic = getticks();
-#ifdef SWIFT_DEBUG_CHECKS
- if (c->hydro.h_max <= 0.) error("Invalid h_max for does_los_intersect");
-#endif
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const struct los_props *LOS_params = e->los_properties;
+ const int verbose = e->verbose;
- /* Distance from LOS to left and bottom cell edge. */
- const double cx_min = c->loc[los->xaxis];
- const double cy_min = c->loc[los->yaxis];
+ /* Start by generating the random sightline positions. */
+ struct line_of_sight *LOS_list = (struct line_of_sight *)malloc(
+ LOS_params->num_tot * sizeof(struct line_of_sight));
- /* Distance from LOS to right and top cell edge. */
- const double cx_max = c->loc[los->xaxis] + c->width[los->xaxis];
- const double cy_max = c->loc[los->yaxis] + c->width[los->yaxis];
+ if (e->nodeID == 0) {
+ generate_sightlines(LOS_list, LOS_params, periodic, dim);
+ if (verbose)
+ message("Generated %i random sightlines.", LOS_params->num_tot);
+ }
- /* Maximum smoothing length of a part in this top level cell. */
- const double hsml = c->hydro.h_max * kernel_gamma;
- const double hsml2 = hsml * hsml;
-
- double dx, dy;
-
- if (los->periodic) {
- dx = min(fabs(nearest(cx_min - los->Xpos, los->dim[los->xaxis])),
- fabs(nearest(cx_max - los->Xpos, los->dim[los->xaxis])));
- dy = min(fabs(nearest(cy_min - los->Ypos, los->dim[los->yaxis])),
- fabs(nearest(cy_max - los->Ypos, los->dim[los->yaxis])));
- } else {
- dx = fabs(cx_max - los->Xpos);
- dy = fabs(cy_max - los->Ypos);
- }
-
- /* Is sightline directly within this top level cell? */
- if (dx < (1.01 * c->width[los->xaxis]) / 2. &&
- dy < (1.01 * c->width[los->yaxis]) / 2.) {
- return 1;
- /* Could a part from this top level cell smooth into the sightline? */
- } else if (dx * dx + dy * dy < hsml2) {
- return 1;
- /* Don't need to work with this top level cell. */
- } else {
- return 0;
- }
-}
-
-/**
- * @brief Main work function for computing line of sights.
- *
- * 1) Construct N random line of sight positions.
- * 2) Loop over each line of sight.
- * - 2.1) Find which top level cells sightline intersects.
- * - 2.2) Loop over each part in these top level cells to see which intersect
- * sightline.
- * - 2.3) Use this count to construct a LOS parts/xparts array.
- * - 2.4) Loop over each part and extract those in sightline to new array.
- * - 2.5) Save sightline parts to HDF5 file.
- *
- * @param e The engine.
- */
-void do_line_of_sight(struct engine *e) {
-
- /* Start counting. */
- const ticks tic = getticks();
-
- const struct space *s = e->s;
- const int periodic = s->periodic;
- const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
- const struct los_props *LOS_params = e->los_properties;
- const int verbose = e->verbose;
-
- /* Start by generating the random sightline positions. */
- struct line_of_sight *LOS_list = (struct line_of_sight *)malloc(
- LOS_params->num_tot * sizeof(struct line_of_sight));
- if (e->nodeID == 0) {
- generate_sightlines(LOS_list, LOS_params, periodic, dim);
- if (verbose)
- message("Generated %i random sightlines.", LOS_params->num_tot);
- }
#ifdef WITH_MPI
+ /* Share the list of LoS with all the MPI ranks */
MPI_Bcast(LOS_list, LOS_params->num_tot * sizeof(struct line_of_sight),
MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
@@ -409,6 +710,7 @@ void do_line_of_sight(struct engine *e) {
h_file = H5Fcreate(fileName, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if (h_file < 0) error("Error while opening file '%s'.", fileName);
}
+
#ifdef WITH_MPI
MPI_Barrier(MPI_COMM_WORLD);
#endif
@@ -420,7 +722,7 @@ void do_line_of_sight(struct engine *e) {
/* Main loop over each random LOS. */
/* ------------------------------- */
- /* Get list of local top level cells. */
+ /* Get list of local non-empty top level cells. */
const struct cell *cells = e->s->cells_top;
const int *local_cells_with_particles = e->s->local_cells_with_particles_top;
const int nr_local_cells_with_particles = s->nr_local_cells_with_particles;
@@ -431,8 +733,7 @@ void do_line_of_sight(struct engine *e) {
/* Create empty top level cell list for this LOS */
int *los_cells_top = (int *)swift_malloc(
"tl_cells_los", nr_local_cells_with_particles * sizeof(int));
- for (int n = 0; n < nr_local_cells_with_particles; n++)
- los_cells_top[n] = 0;
+ bzero(los_cells_top, nr_local_cells_with_particles * sizeof(int));
/* Find all top level cells this LOS will intersect. */
find_intersecting_top_level_cells(e, &LOS_list[j], los_cells_top, cells,
@@ -480,7 +781,9 @@ void do_line_of_sight(struct engine *e) {
MPI_Allgather(&LOS_list[j].particles_in_los_local, 1, MPI_INT, counts, 1,
MPI_INT, MPI_COMM_WORLD);
- for (int k = 0, offset_count = 0; k < e->nr_nodes; k++) {
+ int offset_count = 0;
+ for (int k = 0; k < e->nr_nodes; k++) {
+
/* Total parts in this LOS. */
LOS_list[j].particles_in_los_total += counts[k];
@@ -502,6 +805,10 @@ void do_line_of_sight(struct engine *e) {
message("*WARNING* LOS %i is empty", j);
print_los_info(LOS_list, j);
}
+#ifdef WITH_MPI
+ free(counts);
+ free(offsets);
+#endif
swift_free("tl_cells_los", los_cells_top);
continue;
}
@@ -510,6 +817,7 @@ void do_line_of_sight(struct engine *e) {
struct part *LOS_parts = NULL;
struct xpart *LOS_xparts = NULL;
+ /* Rank 0 allocates more space as it will gather all the data for writing */
if (e->nodeID == 0) {
if ((LOS_parts = (struct part *)swift_malloc(
"los_parts_array",
@@ -535,8 +843,7 @@ void do_line_of_sight(struct engine *e) {
}
/* Loop over each part again, pulling out those in LOS. */
- size_t count = 0;
- double dx, dy, dx2, dy2, hsml, hsml2;
+ int count = 0;
for (int n = 0; n < e->s->nr_local_cells_with_particles; ++n) {
@@ -560,30 +867,31 @@ void do_line_of_sight(struct engine *e) {
continue;
/* Distance from this part to LOS along x dim. */
- dx = cell_parts[i].x[LOS_list[j].xaxis] - LOS_list[j].Xpos;
+ double dx = cell_parts[i].x[LOS_list[j].xaxis] - LOS_list[j].Xpos;
/* Periodic wrap. */
if (LOS_list[j].periodic)
dx = nearest(dx, LOS_list[j].dim[LOS_list[j].xaxis]);
/* Square */
- dx2 = dx * dx;
+ const double dx2 = dx * dx;
/* Smoothing length of this part. */
- hsml = cell_parts[i].h * kernel_gamma;
- hsml2 = hsml * hsml;
+ const double hsml = cell_parts[i].h * kernel_gamma;
+ const double hsml2 = hsml * hsml;
/* Does this part fall into our LOS? */
if (dx2 < hsml2) {
+
/* Distance from this part to LOS along y dim. */
- dy = cell_parts[i].x[LOS_list[j].yaxis] - LOS_list[j].Ypos;
+ double dy = cell_parts[i].x[LOS_list[j].yaxis] - LOS_list[j].Ypos;
/* Periodic wrap. */
if (LOS_list[j].periodic)
dy = nearest(dy, LOS_list[j].dim[LOS_list[j].yaxis]);
/* Square */
- dy2 = dy * dy;
+ const double dy2 = dy * dy;
/* Does this part still fall into our LOS? */
if (dy2 < hsml2) {
@@ -649,7 +957,11 @@ void do_line_of_sight(struct engine *e) {
H5Gclose(h_grp);
}
- /* Free up some memory */
+ /* Free up some memory */
+#ifdef WITH_MPI
+ free(counts);
+ free(offsets);
+#endif
swift_free("tl_cells_los", los_cells_top);
swift_free("los_parts_array", LOS_parts);
swift_free("los_xparts_array", LOS_xparts);
@@ -673,294 +985,6 @@ void do_line_of_sight(struct engine *e) {
clocks_getunit());
}
-/**
- * @brief Write parts in LOS to HDF5 file.
- *
- * @param grp HDF5 group of this LOS.
- * @param j LOS ID.
- * @param N number of parts in this line of sight.
- * @param parts the list of parts in this LOS.
- * @param e The engine.
- * @param xparts the list of xparts in this LOS.
- */
-void write_los_hdf5_datasets(hid_t grp, int j, size_t N,
- const struct part *parts, struct engine *e,
- const struct xpart *xparts) {
-
- /* What kind of run are we working with? */
- struct swift_params *params = e->parameter_file;
- const int with_cosmology = e->policy & engine_policy_cosmology;
- const int with_cooling = e->policy & engine_policy_cooling;
- const int with_temperature = e->policy & engine_policy_temperature;
- const int with_fof = e->policy & engine_policy_fof;
-#ifdef HAVE_VELOCIRAPTOR
- const int with_stf = (e->policy & engine_policy_structure_finding) &&
- (e->s->gpart_group_data != NULL);
-#else
- const int with_stf = 0;
-#endif
-
- int num_fields = 0;
- struct io_props list[100];
-
- /* Find all the gas output fields */
- hydro_write_particles(parts, xparts, list, &num_fields);
- num_fields += chemistry_write_particles(parts, list + num_fields);
- if (with_cooling || with_temperature) {
- num_fields += cooling_write_particles(parts, xparts, list + num_fields,
- e->cooling_func);
- }
- if (with_fof) {
- num_fields += fof_write_parts(parts, xparts, list + num_fields);
- }
- if (with_stf) {
- num_fields += velociraptor_write_parts(parts, xparts, list + num_fields);
- }
- num_fields +=
- tracers_write_particles(parts, xparts, list + num_fields, with_cosmology);
- num_fields +=
- star_formation_write_particles(parts, xparts, list + num_fields);
-
- /* Loop over each output field */
- for (int i = 0; i < num_fields; i++) {
-
- /* Did the user cancel this field? */
- char field[PARSER_MAX_LINE_SIZE];
- sprintf(field, "SelectOutputLOS:%.*s", FIELD_BUFFER_SIZE, list[i].name);
- int should_write = parser_get_opt_param_int(params, field, 1);
-
- /* Write (if selected) */
- if (should_write) write_los_hdf5_dataset(list[i], N, j, e, grp);
- // if (should_write) write_array_single(e, grp, NULL,
- // NULL, NULL,
- // list[i], N,
- // e->internal_units,
- // e->snapshot_units);
- }
-}
-
-/**
- * @brief Writes dataset for a given part attribute.
- *
- * @param p io_props dataset for this attribute.
- * @param N number of parts in this line of sight.
- * @param j Line of sight ID.
- * @param e The engine.
- * @param grp HDF5 group to write to.
- */
-void write_los_hdf5_dataset(const struct io_props props, size_t N, int j,
- struct engine *e, hid_t grp) {
-
- /* Create data space */
- const hid_t h_space = H5Screate(H5S_SIMPLE);
- if (h_space < 0)
- error("Error while creating data space for field '%s'.", props.name);
-
- int rank = 0;
- hsize_t shape[2];
- hsize_t chunk_shape[2];
- if (props.dimension > 1) {
- rank = 2;
- shape[0] = N;
- shape[1] = props.dimension;
- chunk_shape[0] = 1 << 20; /* Just a guess...*/
- chunk_shape[1] = props.dimension;
- } else {
- rank = 1;
- shape[0] = N;
- shape[1] = 0;
- chunk_shape[0] = 1 << 20; /* Just a guess...*/
- chunk_shape[1] = 0;
- }
-
- /* Make sure the chunks are not larger than the dataset */
- if (chunk_shape[0] > N) chunk_shape[0] = N;
-
- /* Change shape of data space */
- hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, shape);
- if (h_err < 0)
- error("Error while changing data space shape for field '%s'.", props.name);
-
- /* Dataset properties */
- const hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
-
- /* Set chunk size */
- h_err = H5Pset_chunk(h_prop, rank, chunk_shape);
- if (h_err < 0)
- error("Error while setting chunk size (%llu, %llu) for field '%s'.",
- chunk_shape[0], chunk_shape[1], props.name);
-
- /* Impose check-sum to verify data corruption */
- h_err = H5Pset_fletcher32(h_prop);
- if (h_err < 0)
- error("Error while setting checksum options for field '%s'.", props.name);
-
- /* Impose data compression */
- if (e->snapshot_compression > 0) {
- h_err = H5Pset_shuffle(h_prop);
- if (h_err < 0)
- error("Error while setting shuffling options for field '%s'.",
- props.name);
-
- h_err = H5Pset_deflate(h_prop, e->snapshot_compression);
- if (h_err < 0)
- error("Error while setting compression options for field '%s'.",
- props.name);
- }
-
- /* Allocate temporary buffer */
- const size_t num_elements = N * props.dimension;
- const size_t typeSize = io_sizeof_type(props.type);
- void *temp = NULL;
- if (swift_memalign("writebuff", (void **)&temp, IO_BUFFER_ALIGNMENT,
- num_elements * typeSize) != 0)
- error("Unable to allocate temporary i/o buffer");
-
- /* Copy particle data to temp buffer */
- io_copy_temp_buffer(temp, e, props, N, e->internal_units, e->snapshot_units);
-
- /* Create dataset */
- char att_name[200];
- sprintf(att_name, "/LOS_%04i/%s", j, props.name);
- const hid_t h_data = H5Dcreate(grp, att_name, io_hdf5_type(props.type),
- h_space, H5P_DEFAULT, h_prop, H5P_DEFAULT);
- if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
-
- /* Write dataset */
- herr_t status = H5Dwrite(h_data, io_hdf5_type(props.type), H5S_ALL, H5S_ALL,
- H5P_DEFAULT, temp);
- if (status < 0) error("Error while writing data array '%s'.", props.name);
-
- /* Write unit conversion factors for this data set */
- char buffer[FIELD_BUFFER_SIZE] = {0};
- units_cgs_conversion_string(buffer, e->snapshot_units, props.units,
- props.scale_factor_exponent);
- float baseUnitsExp[5];
- units_get_base_unit_exponents_array(baseUnitsExp, props.units);
- io_write_attribute_f(h_data, "U_M exponent", baseUnitsExp[UNIT_MASS]);
- io_write_attribute_f(h_data, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
- io_write_attribute_f(h_data, "U_t exponent", baseUnitsExp[UNIT_TIME]);
- io_write_attribute_f(h_data, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
- io_write_attribute_f(h_data, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
- io_write_attribute_f(h_data, "h-scale exponent", 0.f);
- io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
- io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
-
- /* Write the actual number this conversion factor corresponds to */
- const double factor =
- units_cgs_conversion_factor(e->snapshot_units, props.units);
- io_write_attribute_d(
- h_data,
- "Conversion factor to CGS (not including cosmological corrections)",
- factor);
- io_write_attribute_d(
- h_data,
- "Conversion factor to physical CGS (including cosmological corrections)",
- factor * pow(e->cosmology->a, props.scale_factor_exponent));
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (strlen(props.description) == 0)
- error("Invalid (empty) description of the field '%s'", props.name);
-#endif
-
- /* Write the full description */
- io_write_attribute_s(h_data, "Description", props.description);
-
- /* Free and close everything */
- swift_free("writebuff", temp);
- H5Pclose(h_prop);
- H5Dclose(h_data);
- H5Sclose(h_space);
-}
-
-/**
- * @brief Writes HDF5 headers and information groups for this line of sight.
- *
- * @param h_file HDF5 file reference.
- * @param e The engine.
- * @param LOS_params The line of sight params.
- */
-void write_hdf5_header(hid_t h_file, const struct engine *e,
- const struct los_props *LOS_params,
- const size_t total_num_parts_in_los) {
- /* Open header to write simulation properties */
- hid_t h_grp =
- H5Gcreate(h_file, "/Header", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating file header\n");
-
- /* Convert basic output information to snapshot units */
- const double factor_time = units_conversion_factor(
- e->internal_units, e->snapshot_units, UNIT_CONV_TIME);
- const double factor_length = units_conversion_factor(
- e->internal_units, e->snapshot_units, UNIT_CONV_LENGTH);
- const double dblTime = e->time * factor_time;
- const double dim[3] = {e->s->dim[0] * factor_length,
- e->s->dim[1] * factor_length,
- e->s->dim[2] * factor_length};
-
- io_write_attribute(h_grp, "BoxSize", DOUBLE, dim, 3);
- io_write_attribute(h_grp, "Time", DOUBLE, &dblTime, 1);
- const int dimension = (int)hydro_dimension;
- io_write_attribute(h_grp, "Dimension", INT, &dimension, 1);
- io_write_attribute(h_grp, "Redshift", DOUBLE, &e->cosmology->z, 1);
- io_write_attribute(h_grp, "Scale-factor", DOUBLE, &e->cosmology->a, 1);
- io_write_attribute_s(h_grp, "Code", "SWIFT");
- io_write_attribute_s(h_grp, "RunName", e->run_name);
- const int num_files_per_snapshot = 1;
- io_write_attribute(h_grp, "NumFilesPerSnapshot", INT, &num_files_per_snapshot,
- 1);
- io_write_attribute_s(h_grp, "OutputType", "LineOfSight");
-
- /* GADGET-2 legacy values */
- /* Number of particles of each type */
- long long N_total[swift_type_count] = {0};
- N_total[0] = total_num_parts_in_los;
- unsigned int numParticles[swift_type_count] = {0};
- unsigned int numParticlesHighWord[swift_type_count] = {0};
- for (int ptype = 0; ptype < swift_type_count; ++ptype) {
- numParticles[ptype] = (unsigned int)N_total[ptype];
- numParticlesHighWord[ptype] = (unsigned int)(N_total[ptype] >> 32);
- }
- io_write_attribute(h_grp, "NumPart_ThisFile", LONGLONG, N_total,
- swift_type_count);
- io_write_attribute(h_grp, "NumPart_Total", UINT, numParticles,
- swift_type_count);
- io_write_attribute(h_grp, "NumPart_Total_HighWord", UINT,
- numParticlesHighWord, swift_type_count);
-
- /* Store the time at which the snapshot was written */
- time_t tm = time(NULL);
- struct tm *timeinfo = localtime(&tm);
- char snapshot_date[64];
- strftime(snapshot_date, 64, "%T %F %Z", timeinfo);
- io_write_attribute_s(h_grp, "Snapshot date", snapshot_date);
-
- /* Close group */
- H5Gclose(h_grp);
-
- io_write_meta_data(h_file, e, e->internal_units, e->snapshot_units);
-
- /* Print the LOS properties */
- h_grp = H5Gcreate(h_file, "/LineOfSightParameters", H5P_DEFAULT, H5P_DEFAULT,
- H5P_DEFAULT);
- if (h_grp < 0) error("Error while creating LOS group");
-
- /* Record this LOS attributes */
- const int num_los_per_axis[3] = {LOS_params->num_along_x,
- LOS_params->num_along_y,
- LOS_params->num_along_z};
- io_write_attribute(h_grp, "NumLineOfSight_PerAxis", INT, num_los_per_axis, 3);
- io_write_attribute(h_grp, "NumLineOfSight_Total", INT, &LOS_params->num_tot,
- 1);
- io_write_attribute(h_grp, "AllowedLOSRangeX", DOUBLE,
- LOS_params->allowed_losrange_x, 2);
- io_write_attribute(h_grp, "AllowedLOSRangeY", DOUBLE,
- LOS_params->allowed_losrange_y, 2);
- io_write_attribute(h_grp, "AllowedLOSRangeZ", DOUBLE,
- LOS_params->allowed_losrange_z, 2);
- H5Gclose(h_grp);
-}
-
/**
* @brief Write a los_props struct to the given FILE as a stream of bytes.
*
diff --git a/src/line_of_sight.h b/src/line_of_sight.h
index 4fe9eb12779eeee289b3b09247c597633c04666e..0e79a9a50126a76cee152d9fc9612b7b33961407 100644
--- a/src/line_of_sight.h
+++ b/src/line_of_sight.h
@@ -26,8 +26,8 @@
#include "engine.h"
#include "io_properties.h"
-/*
- * Maps the LOS axis geometry to the simulation axis geometry.
+/**
+ * @brief Maps the LOS axis geometry to the simulation axis geometry.
*
* Sigtlines will always shoot down the los_direction_z,
* which can map to x,y or z of the simulation geometry.
@@ -36,106 +36,87 @@
* the plane orthogonal to the LOS direction. The random
* sightline positions are created on this plane.
*/
-enum los_direction {
- simulation_x_axis = 0,
- simulation_y_axis = 1,
- simulation_z_axis = 2
-};
+enum los_direction { simulation_x_axis, simulation_y_axis, simulation_z_axis };
+/**
+ * @brief Properties of a single line-of-sight
+ */
struct line_of_sight {
- /* Simulation axis the LOS shoots down. */
+
+ /*! Simulation axis the LOS shoots down. */
enum los_direction zaxis;
- /* The two remaining axes defining the plane orthogonal to the sightline. */
+ /*! The two remaining axes defining the plane orthogonal to the sightline. */
enum los_direction xaxis, yaxis;
- /* Sightline position along los_direction_x. */
+ /*! Sightline position along los_direction_x. */
double Xpos;
- /* Sightline position along los_direction_y. */
+ /*! Sightline position along los_direction_y. */
double Ypos;
- /* Number of parts in LOS. */
- size_t particles_in_los_total;
+ /*! Number of parts in LOS. */
+ int particles_in_los_total;
- /* Number of parts in LOS on this node. */
- size_t particles_in_los_local;
+ /*! Number of parts in LOS on this node. */
+ int particles_in_los_local;
- /* Is the simulation periodic? */
+ /*! Is the simulation periodic? */
int periodic;
- /* Dimensions of the space. */
+ /*! Dimensions of the space. */
double dim[3];
- /* How many top level cells does ths LOS intersect? */
+ /*! How many top level cells does ths LOS intersect? */
int num_intersecting_top_level_cells;
- /* The min--max range to consider for parts in LOS. */
+ /*! The min--max range to consider for parts in LOS. */
double range_when_shooting_down_axis[2];
};
+/**
+ * @brief Properties of the line-of-sight computation
+ */
struct los_props {
- /* Number of sightlines shooting down simulation z axis. */
+
+ /*! Number of sightlines shooting down simulation z axis. */
int num_along_z;
- /* Number of sightlines shooting down simulation x axis. */
+ /*! Number of sightlines shooting down simulation x axis. */
int num_along_x;
- /* Number of sightlines shooting down simulation y axis. */
+ /*! Number of sightlines shooting down simulation y axis. */
int num_along_y;
- /* Total number of sightlines. */
+ /*! Total number of sightlines. */
int num_tot;
- /* The min--max range along the simulation x axis random sightlines are
+ /*! The min--max range along the simulation x axis random sightlines are
* allowed. */
double allowed_losrange_x[2];
- /* The min--max range along the simulation y axis random sightlines are
+ /*! The min--max range along the simulation y axis random sightlines are
* allowed. */
double allowed_losrange_y[2];
- /* The min--max range along the simulation z axis random sightlines are
+ /*! The min--max range along the simulation z axis random sightlines are
* allowed. */
double allowed_losrange_z[2];
- /* The min--max range to consider when LOS is shooting down each
- * simulation axis. */
+ /*! The min--max range to consider when LOS is shooting down each simulation
+ * axis. */
double range_when_shooting_down_axis[3][2];
- /* Basename for line of sight HDF5 files. */
+ /*! Base name for line of sight HDF5 files. */
char basename[200];
};
-double los_periodic(double x, double dim);
-void generate_sightlines(struct line_of_sight *Los,
- const struct los_props *params, const int periodic,
- const double dim[3]);
void print_los_info(const struct line_of_sight *Los, const int i);
void do_line_of_sight(struct engine *e);
-void los_init(double dim[3], struct los_props *los_params,
+void los_init(const double dim[3], struct los_props *los_params,
struct swift_params *params);
-void write_los_hdf5_datasets(hid_t grp, int j, size_t N,
- const struct part *parts, struct engine *e,
- const struct xpart *xparts);
-void write_los_hdf5_dataset(const struct io_props p, size_t N, int j,
- struct engine *e, hid_t grp);
-void write_hdf5_header(hid_t h_file, const struct engine *e,
- const struct los_props *LOS_params,
- const size_t total_num_parts_in_los);
-void create_sightline(const double Xpos, const double Ypos,
- enum los_direction xaxis, enum los_direction yaxis,
- enum los_direction zaxis, const int periodic,
- const double dim[3], struct line_of_sight *los,
- const double range_when_shooting_down_axis[2]);
+
void los_struct_dump(const struct los_props *internal_los, FILE *stream);
void los_struct_restore(const struct los_props *internal_los, FILE *stream);
-int does_los_intersect(const struct cell *c, const struct line_of_sight *los);
-void find_intersecting_top_level_cells(const struct engine *e,
- struct line_of_sight *los,
- int *los_cells_top,
- const struct cell *cells,
- const int *local_cells_with_particles,
- const int nr_local_cells_with_particles);
#endif /* SWIFT_LOS_H */
diff --git a/src/part.c b/src/part.c
index 91234f6215276f51b5a01b0bc9b0b34d7f621484..be685f1d862e5fe69fd45c73998fbfc0a0b3b2ff 100644
--- a/src/part.c
+++ b/src/part.c
@@ -31,7 +31,6 @@
/* Local headers */
#include "error.h"
#include "hydro.h"
-#include "line_of_sight.h"
#include "threadpool.h"
/**
diff --git a/src/single_io.h b/src/single_io.h
index 3c9bfbae637dbc9832de9c71340224a30cbfa3cd..db9ac92b8b1cc71e588391a2ff0a78839e3e61ec 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -29,6 +29,7 @@
struct engine;
struct unit_system;
+struct io_props;
void read_ic_single(const char* fileName,
const struct unit_system* internal_units, double dim[3],
@@ -47,13 +48,12 @@ void write_output_single(struct engine* e,
#endif /* HAVE_HDF5 && !WITH_MPI */
#ifdef HAVE_HDF5
-#include "io_properties.h"
-
void write_array_single(const struct engine* e, hid_t grp, char* fileName,
FILE* xmfFile, char* partTypeGroupName,
const struct io_props props, size_t N,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units);
-#endif
+
+#endif /* HAVE_HDF5 */
#endif /* SWIFT_SINGLE_IO_H */