diff --git a/examples/main.c b/examples/main.c
index b3e85987d753d1f54c01dc68e64efce12d34d4f5..edf9610a2ea97acfb0362a0954b3abe941ca31ef 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -106,6 +106,7 @@ int main(int argc, char *argv[]) {
struct spart *sparts = NULL;
struct bpart *bparts = NULL;
struct unit_system us;
+ struct los_props los_properties;
int nr_nodes = 1, myrank = 0;
@@ -1063,6 +1064,9 @@ int main(int argc, char *argv[]) {
with_hydro, with_self_gravity, with_star_formation,
with_DM_background_particles, talking, dry_run, nr_nodes);
+ /* Initialise the line of sight properties. */
+ los_init(s.dim, &los_properties, params);
+
if (myrank == 0) {
clocks_gettime(&toc);
message("space_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -1196,7 +1200,7 @@ int main(int argc, char *argv[]) {
&reparttype, &us, &prog_const, &cosmo, &hydro_properties,
&entropy_floor, &gravity_properties, &stars_properties,
&black_holes_properties, &feedback_properties, &mesh, &potential,
- &cooling_func, &starform, &chemistry, &fof_properties);
+ &cooling_func, &starform, &chemistry, &fof_properties, &los_properties);
engine_config(/*restart=*/0, /*fof=*/0, &e, params, nr_nodes, myrank,
nr_threads, with_aff, talking, restart_file);
diff --git a/src/engine.c b/src/engine.c
index 1204095786444ae92a22f5e3a6e0ae763018e19b..856872f60ce0ab1a820fc3dbdbd2b757281bd5f3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2805,8 +2805,12 @@ void engine_check_for_dumps(struct engine *e) {
}
/* Do we want to write a line of sight file? */
- ti_output = e->ti_current;
- type = output_los;
+ if (e->ti_end_min > e->ti_next_los && e->ti_next_los > 0) {
+ if (e->ti_next_los < ti_output) {
+ ti_output = e->ti_next_los;
+ type = output_los;
+ }
+ }
/* Store information before attempting extra dump-related drifts */
const integertime_t ti_current = e->ti_current;
@@ -2897,6 +2901,8 @@ void engine_check_for_dumps(struct engine *e) {
case output_los:
do_line_of_sight(e);
+ engine_compute_next_los_time(e);
+
break;
default:
@@ -2935,6 +2941,14 @@ void engine_check_for_dumps(struct engine *e) {
}
}
+ /* Do line of sight ? */
+ if (e->ti_end_min > e->ti_next_los && e->ti_next_los > 0) {
+ if (e->ti_next_los < ti_output) {
+ ti_output = e->ti_next_los;
+ type = output_los;
+ }
+ }
+
} /* While loop over output types */
/* Restore the information we stored */
@@ -3835,7 +3849,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
- struct fof_props *fof_properties) {
+ struct fof_props *fof_properties, struct los_props *los_properties) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -3884,6 +3898,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
units_init_default(e->snapshot_units, params, "Snapshots", internal_units);
e->snapshot_output_count = 0;
e->stf_output_count = 0;
+ e->los_output_count = 0;
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->dt_max_RMS_displacement = FLT_MAX;
@@ -3916,6 +3931,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->fof_properties = fof_properties;
e->parameter_file = params;
e->stf_this_timestep = 0;
+ e->los_properties = los_properties;
#ifdef WITH_MPI
e->usertime_last_step = 0.0;
e->systime_last_step = 0.0;
@@ -4472,6 +4488,9 @@ void engine_config(int restart, int fof, struct engine *e,
/* Find the time of the first statistics output */
engine_compute_next_statistics_time(e);
+ /* Find the time of the first line of sight output */
+ engine_compute_next_los_time(e);
+
/* Find the time of the first stf output */
if (e->policy & engine_policy_structure_finding) {
engine_compute_next_stf_time(e);
@@ -4884,6 +4903,77 @@ void engine_compute_next_statistics_time(struct engine *e) {
}
}
+/**
+ * @brief Computes the next time (on the time line) for a line of sight dump
+ *
+ * @param e The #engine.
+ */
+void engine_compute_next_los_time(struct engine *e) {
+ /* Do output_list file case */
+ if (e->output_list_los) {
+ output_list_read_next_time(e->output_list_los, e, "line of sights",
+ &e->ti_next_los);
+ return;
+ }
+
+ /* Find upper-bound on last output */
+ double time_end;
+ if (e->policy & engine_policy_cosmology)
+ time_end = e->cosmology->a_end * e->delta_time_los;
+ else
+ time_end = e->time_end + e->delta_time_los;
+
+ /* Find next los above current time */
+ double time;
+ if (e->policy & engine_policy_cosmology)
+ time = e->a_first_los;
+ else
+ time = e->time_first_los;
+
+ int found_los_time = 0;
+ while (time < time_end) {
+
+ /* Output time on the integer timeline */
+ if (e->policy & engine_policy_cosmology)
+ e->ti_next_los = log(time / e->cosmology->a_begin) / e->time_base;
+ else
+ e->ti_next_los = (time - e->time_begin) / e->time_base;
+
+ /* Found it? */
+ if (e->ti_next_los > e->ti_current) {
+ found_los_time = 1;
+ break;
+ }
+
+ if (e->policy & engine_policy_cosmology)
+ time *= e->delta_time_los;
+ else
+ time += e->delta_time_los;
+ }
+
+ /* Deal with last line of sight */
+ if (!found_los_time) {
+ e->ti_next_los = -1;
+ if (e->verbose) message("No further LOS output time.");
+ } else {
+
+ /* Be nice, talk... */
+ if (e->policy & engine_policy_cosmology) {
+ const double next_los_time =
+ exp(e->ti_next_los * e->time_base) * e->cosmology->a_begin;
+ if (e->verbose)
+ message("Next output time for line of sight set to a=%e.",
+ next_los_time);
+ } else {
+ const double next_los_time =
+ e->ti_next_los * e->time_base + e->time_begin;
+ if (e->verbose)
+ message("Next output time for line of sight set to t=%e.",
+ next_los_time);
+ }
+ }
+}
+
/**
* @brief Computes the next time (on the time line) for structure finding
*
@@ -5057,6 +5147,19 @@ void engine_init_output_lists(struct engine *e, struct swift_params *params) {
else
e->time_first_stf_output = stf_time_first;
}
+
+ /* Deal with line of sight */
+ double los_time_first;
+ e->output_list_los = NULL;
+ output_list_init(&e->output_list_los, e, "LineOfSight",
+ &e->delta_time_los, &los_time_first);
+
+ if (e->output_list_los) {
+ if (e->policy & engine_policy_cosmology)
+ e->a_first_los = stf_time_first;
+ else
+ e->time_first_los = stf_time_first;
+ }
}
/**
diff --git a/src/engine.h b/src/engine.h
index 3543667a563587e126cf191968d93c7c5e7b4e0a..e8d69f8e060f5b38fed3ffafe4135e0af51e6611 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -488,6 +488,17 @@ struct engine {
/* Has there been an stf this timestep? */
char stf_this_timestep;
+ /* Line of sight properties. */
+ struct los_props* los_properties;
+
+ /* Line of sight outputs information. */
+ struct output_list *output_list_los;
+ double a_first_los;
+ double time_first_los;
+ double delta_time_los;
+ integertime_t ti_next_los;
+ int los_output_count;
+
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Run brute force checks only on steps when all gparts active? */
int force_checks_only_all_active;
@@ -510,6 +521,7 @@ void engine_compute_next_snapshot_time(struct engine *e);
void engine_compute_next_stf_time(struct engine *e);
void engine_compute_next_fof_time(struct engine *e);
void engine_compute_next_statistics_time(struct engine *e);
+void engine_compute_next_los_time(struct engine *e);
void engine_recompute_displacement_constraint(struct engine *e);
void engine_unskip(struct engine *e);
void engine_unskip_timestep_communications(struct engine *e);
@@ -538,7 +550,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
struct cooling_function_data *cooling_func,
const struct star_formation *starform,
const struct chemistry_global_data *chemistry,
- struct fof_props *fof_properties);
+ struct fof_props *fof_properties,
+ struct los_props *los_properties);
void engine_config(int restart, int fof, struct engine *e,
struct swift_params *params, int nr_nodes, int nodeID,
int nr_threads, int with_aff, int verbose,
diff --git a/src/line_of_sight.c b/src/line_of_sight.c
index 5fa1981b677e4a26fc599ed501a9f3b398960d49..3df5e1002721d0969db64c13a20f8d0244f8d571 100644
--- a/src/line_of_sight.c
+++ b/src/line_of_sight.c
@@ -14,6 +14,36 @@
#include "line_of_sight.h"
#include "periodic.h"
+void los_init(double dim[3], struct los_props *los_params,
+ struct swift_params *params) {
+ /* How many line of sights in each plane. */
+ los_params->num_along_xy =
+ parser_get_opt_param_int(params, "LineOfSight:num_along_xy", 0);
+ los_params->num_along_yz =
+ parser_get_opt_param_int(params, "LineOfSight:num_along_yz", 0);
+ los_params->num_along_xz =
+ parser_get_opt_param_int(params, "LineOfSight:num_along_xz", 0);
+
+ /* Min/max range across x,y and z where random LOS's are allowed. */
+ los_params->xmin =
+ parser_get_opt_param_double(params, "LineOfSight:xmin", 0.);
+ los_params->xmax =
+ parser_get_opt_param_double(params, "LineOfSight:xmax", dim[0]);
+ los_params->ymin =
+ parser_get_opt_param_double(params, "LineOfSight:ymin", 0.);
+ los_params->ymax =
+ parser_get_opt_param_double(params, "LineOfSight:ymax", dim[1]);
+ los_params->zmin =
+ parser_get_opt_param_double(params, "LineOfSight:zmin", 0.);
+ los_params->zmax =
+ parser_get_opt_param_double(params, "LineOfSight:zmax", dim[2]);
+
+ /* Compute total number of sightlines. */
+ los_params->num_tot = los_params->num_along_xy +
+ los_params->num_along_yz +
+ los_params->num_along_xz;
+}
+
/**
* @brief Generates random sightline positions.
*
@@ -23,7 +53,7 @@
* @param params Sightline parameters.
*/
void generate_line_of_sights(struct line_of_sight *Los,
- struct line_of_sight_params *params) {
+ struct los_props *params) {
/* Keep track of number of sightlines. */
int count = 0;
@@ -84,7 +114,7 @@ void generate_line_of_sights(struct line_of_sight *Los,
* @param params Sightline parameters.
*/
void print_los_info(struct line_of_sight *Los,
- struct line_of_sight_params *params) {
+ struct los_props *params) {
printf("\nPrinting LOS information...\n");
for (int i = 0; i < params->num_tot; i++) {
@@ -100,26 +130,14 @@ void do_line_of_sight(struct engine *e) {
const struct space *s = e->s;
const size_t nr_parts = s->nr_parts;
-
- /* LOS params, dummy until included in parameter file. */
- struct line_of_sight_params LOS_params = {.num_along_xy = 1,
- .num_along_yz = 1,
- .num_along_xz = 1,
- .xmin = 0,
- .xmax = 10,
- .ymin = 0,
- .ymax = 10,
- .zmin = 0,
- .zmax = 10};
- LOS_params.num_tot = LOS_params.num_along_xy + LOS_params.num_along_yz +
- LOS_params.num_along_xz;
+ struct los_props *LOS_params = e->los_properties;
/* 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_line_of_sights(LOS_list, &LOS_params);
+ LOS_params->num_tot * sizeof(struct line_of_sight));
+ if (e->nodeID == 0) generate_line_of_sights(LOS_list, LOS_params);
#ifdef WITH_MPI
- MPI_Bcast(LOS_list, LOS_params.num_tot * sizeof(struct line_of_sight),
+ MPI_Bcast(LOS_list, LOS_params->num_tot * sizeof(struct line_of_sight),
MPI_BYTE, 0, MPI_COMM_WORLD);
#endif
@@ -129,7 +147,7 @@ void do_line_of_sight(struct engine *e) {
if (e->nodeID == 0) {
char filename[200];
- sprintf(filename, "los_%04i.csv", engine_current_step);
+ sprintf(filename, "los_%04i.csv", e->los_output_count);
f = fopen(filename, "w");
if (f == NULL) error("Error opening los file.");
@@ -138,7 +156,7 @@ void do_line_of_sight(struct engine *e) {
}
/* Loop over each random LOS. */
- for (int j = 0; j < LOS_params.num_tot; j++) {
+ for (int j = 0; j < LOS_params->num_tot; j++) {
/* Loop over each gas particle to find those in LOS. */
for (size_t i = 0; i < nr_parts; i++) {
@@ -277,9 +295,12 @@ void do_line_of_sight(struct engine *e) {
if (e->nodeID == 0) fclose(f);
//#ifdef SWIFT_DEBUG_CHECKS
- if (e->nodeID == 0) print_los_info(LOS_list, &LOS_params);
+ if (e->nodeID == 0) print_los_info(LOS_list, LOS_params);
//#endif
-
+
+ /* Up the count. */
+ e->los_output_count++;
+
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
}
diff --git a/src/line_of_sight.h b/src/line_of_sight.h
index e349e96b5512fcf0c13a1113b599dcaccc95bda7..f7d92644ad00c694856d9b30f5abadd856bf0c8a 100644
--- a/src/line_of_sight.h
+++ b/src/line_of_sight.h
@@ -13,7 +13,7 @@ struct line_of_sight_particles {
float h;
};
-struct line_of_sight_params {
+struct los_props {
int num_along_xy;
int num_along_yz;
int num_along_xz;
@@ -26,11 +26,10 @@ struct line_of_sight_params {
};
double los_periodic(double x, double dim);
-
void generate_line_of_sights(struct line_of_sight *Los,
- struct line_of_sight_params *params);
-
+ struct los_props *params);
void print_los_info(struct line_of_sight *Los,
- struct line_of_sight_params *params);
-
+ struct los_props *params);
void do_line_of_sight(struct engine *e);
+void los_init(double dim[3], struct los_props *los_params,
+ struct swift_params *params);
diff --git a/src/swift.h b/src/swift.h
index c7a5f8d5acbcc17dd467cbcfb914eac025894dae..3c659aa7dce1d78d9887187da31e4ee500ad6848 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -50,6 +50,7 @@
#include "hashmap.h"
#include "hydro.h"
#include "hydro_properties.h"
+#include "line_of_sight.h"
#include "lock.h"
#include "logger.h"
#include "logger_io.h"