diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml index 92323a77ce08d7c3b67ec0bd243ae4bec263af16..093f82ada721a0ae2854acce8ebc7719ed4c792b 100644 --- a/examples/EAGLE_6/eagle_6.yml +++ b/examples/EAGLE_6/eagle_6.yml @@ -9,10 +9,9 @@ InternalUnitSystem: # Structure finding on the fly options StructureFinding: config_file_name: stf_input.cfg - output_file_name: stf_output.out - output_time_format: 0 - output_times: 100 - + output_file_name: ./halo/stf + output_time_format: 1 + output_times: 1e-5 # Define the system of units to use int VELOCIraptor. VelociraptorUnitSystem: diff --git a/src/engine.c b/src/engine.c index 36be6ddf89aa7a866f8f557f96c385298505a64b..25223332c0c8d7e07b3fd6e956eb107bbd767fdc 100644 --- a/src/engine.c +++ b/src/engine.c @@ -4489,6 +4489,15 @@ void engine_step(struct engine *e) { if (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) e->dump_snapshot = 1; + /* Do we want to perform structure finding? */ + if ((e->policy & engine_policy_structure_finding)) { + double time = e->timeFirstSnapshot; + if(e->stf_output_time_format == 0 && e->step%(int)e->delta_time_stf_freq == 0) + e->run_stf = 1; + else if(e->stf_output_time_format == 1 && e->ti_end_min >= e->ti_nextSTF && e->ti_nextSTF > 0) + e->run_stf = 1; + } + /* Drift everybody (i.e. what has not yet been drifted) */ /* to the current time */ int drifted_all = @@ -4521,8 +4530,15 @@ void engine_step(struct engine *e) { e->step_props |= engine_step_prop_statistics; } - /* Invoke VELOCIraptor every 250 timesteps thereafter. */ - if ((e->policy & engine_policy_structure_finding) && e->step%250 == 0) velociraptor_invoke(e); + /* Perform structure finding? */ + if (e->run_stf) { + velociraptor_invoke(e); + + /* ... and find the next output time */ + engine_compute_next_stf_time(e); + + e->run_stf = 0; + } /* Now apply all the collected time step updates and particle counts. */ collectgroup1_apply(&e->collect_group1, e); @@ -5341,11 +5357,25 @@ void engine_config(int restart, struct engine *e, e->restart_file = restart_file; e->restart_next = 0; e->restart_dt = 0; + e->delta_time_stf_freq = 0; + e->stf_output_time_format = 0; + e->ti_nextSTF = 0; + e->run_stf = 0; engine_rank = nodeID; /* Initialise VELOCIraptor. */ - if (e->policy & engine_policy_structure_finding) velociraptor_init(e); - + if (e->policy & engine_policy_structure_finding) { + velociraptor_init(e); + e->stf_output_time_format = parser_get_param_int(params, "StructureFinding:output_time_format"); + if(e->stf_output_time_format == 0) { + e->delta_time_stf_freq = (double)parser_get_param_int(params, "StructureFinding:output_times"); + } + else if(e->stf_output_time_format == 1) { + e->delta_time_stf_freq = parser_get_param_double(params, "StructureFinding:output_times"); + } + else error("Invalid flag (%d) set for output time format of structure finding.", e->stf_output_time_format); + } + /* Get the number of queues */ int nr_queues = parser_get_opt_param_int(params, "Scheduler:nr_queues", nr_threads); @@ -5612,8 +5642,13 @@ void engine_config(int restart, struct engine *e, "Time of first snapshot (%e) must be after the simulation start t=%e.", e->timeFirstSnapshot, e->time_begin); + /* Find the time of the first stf output */ + //engine_compute_next_stf_time(e); + message("Next STF step will be: %d", e->ti_nextSTF); + /* Find the time of the first output */ engine_compute_next_snapshot_time(e); + message("Next snapshot step will be: %d", e->ti_nextSnapshot); /* Whether restarts are enabled. Yes by default. Can be changed on restart. */ e->restart_dump = parser_get_opt_param_int(params, "Restarts:enable", 1); @@ -5843,6 +5878,65 @@ void engine_compute_next_snapshot_time(struct engine *e) { } } +/** + * @brief Computes the next time (on the time line) for structure finding + * + * @param e The #engine. + */ +void engine_compute_next_stf_time(struct engine *e) { + + message("dt: %lf", e->delta_time_stf_freq); + /* Find upper-bound on last output */ + double time_end; + if (e->policy & engine_policy_cosmology) + time_end = e->cosmology->a_end * e->delta_time_stf_freq; + else + time_end = e->time_end + e->delta_time_stf_freq; + + /* Find next snasphot above current time */ + double time = e->timeFirstSnapshot; + + message("time: %lf, time_end: %lf", time, e->time_end); + message("dt: %lf", e->delta_time_stf_freq); + + while (time < time_end) { + + /* Output time on the integer timeline */ + if (e->policy & engine_policy_cosmology) + e->ti_nextSTF = log(time / e->cosmology->a_begin) / e->time_base; + else + e->ti_nextSTF = (time - e->time_begin) / e->time_base; + + /* Found it? */ + if (e->ti_nextSTF > e->ti_current) break; + + if (e->policy & engine_policy_cosmology) + time *= e->delta_time_stf_freq; + else + time += e->delta_time_stf_freq; + } + + /* Deal with last snapshot */ + if (e->ti_nextSTF >= max_nr_timesteps) { + e->ti_nextSTF = -1; + if (e->verbose) message("No further output time."); + } else { + + /* Be nice, talk... */ + if (e->policy & engine_policy_cosmology) { + const float next_snapshot_time = + exp(e->ti_nextSTF * e->time_base) * e->cosmology->a_begin; + if (e->verbose) + message("Next output time set to a=%e.", next_snapshot_time); + } else { + const float next_snapshot_time = + e->ti_nextSTF * e->time_base + e->time_begin; + if (e->verbose) + message("Next output time set to t=%e.", next_snapshot_time); + } + } +} + /** * @brief Frees up the memory allocated for this #engine */ diff --git a/src/engine.h b/src/engine.h index f74659ae9e56c3c52a7c0e5dbee04d7b0b9c0106..5da122ce96c99c918fd64e4b13348418a402dfcd 100644 --- a/src/engine.h +++ b/src/engine.h @@ -201,6 +201,12 @@ struct engine { struct unit_system *snapshotUnits; int snapshotOutputCount; + /* Structure finding information */ + int run_stf; + int stf_output_time_format; + double delta_time_stf_freq; + integertime_t ti_nextSTF; + /* Statistics information */ FILE *file_stats; double timeLastStatistics; @@ -325,6 +331,7 @@ struct engine { /* Function prototypes. */ void engine_barrier(struct engine *e); void engine_compute_next_snapshot_time(struct engine *e); +void engine_compute_next_stf_time(struct engine *e); void engine_unskip(struct engine *e); void engine_drift_all(struct engine *e); void engine_drift_top_multipoles(struct engine *e);