diff --git a/doc/RTD/source/ParameterFiles/parameter_description.rst b/doc/RTD/source/ParameterFiles/parameter_description.rst index 9af0ba21a5c318eb1126f0e9d7f6e19208f80b94..aca912da99771085d9d1068adf2f1587980c013d 100644 --- a/doc/RTD/source/ParameterFiles/parameter_description.rst +++ b/doc/RTD/source/ParameterFiles/parameter_description.rst @@ -915,6 +915,21 @@ output as a zero-padded integer. For example, if the base-name is "eagle" and snapshot 7 was just dumped, with ``dump_command`` set to ``./postprocess.sh``, then SWIFT will run ``./postprocess.sh eagle 0007``. +For some quantities, especially in the subgrid models, it can be advantageous to +start recording numbers at a fixed time before the dump of a snapshot. Classic +examples are an averaged star-formation rate or accretion rate onto BHs. For the +subgrid models that support it, the triggers can be specified by setting the +following parameters: + +* for gas: ``recording_triggers_part`` (no default, array of size set by each subgrid model) +* for stars: ``recording_triggers_spart`` (no default, array of size set by each subgrid model) +* for BHs: ``recording_triggers_bpart`` (no default, array of size set by each subgrid model) + +The time is specified in internal time units (See the :ref:`Parameters_units` +section) and a recording can be ignored by setting the parameter to ``-1``. Note +that the code will verify that the recording time is smaller than the gap in +between consecutive snapshot dumps. + Finally, it is possible to specify a different system of units for the snapshots than the one that was used internally by SWIFT. The format is identical to the one described above (See the :ref:`Parameters_units` section) and read: @@ -954,7 +969,7 @@ would have: invoke_fof: 1 compression: 3 distributed: 1 - lustre_OST_count: 48 # System has 48 Lustre OSTs to distribute the files over + lustre_OST_count: 48 # System has 48 Lustre OSTs to distribute the files over UnitLength_in_cgs: 1. # Use cm in outputs UnitMass_in_cgs: 1. # Use grams in outputs UnitVelocity_in_cgs: 1. # Use cm/s in outputs @@ -964,6 +979,12 @@ would have: subsample_fraction: [0, 0.01, 0, 0, 0, 0, 0.1] # Write 1% of the DM parts and 10% of the neutrinos run_on_dump: 1 dump_command: ./submit_analysis.sh + use_delta_from_edge: 1 + delta_from_edge: 1e-6 # Move particles away from the edge by 1-e6 of the length unit. + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_spart: [-1, -1] # No recording + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + Some additional specific options for the snapshot outputs are described in the following pages: diff --git a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml index 118cf8e1556217dbc249a4e3ffa9665d01354cb9..55d97e3376553d880d534f0be7bb96b1080e17dd 100644 --- a/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml +++ b/examples/EAGLE_ICs/EAGLE_100/eagle_100.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml index 89c67559649181485e8c52bd8e0f048e3bdb146d..d37e62d290d63712aad0daa250e0e4289045beee 100644 --- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml +++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml index f8e0bf7042b28a783c2179d44e3fd4089cc6ad42..739f4a2a769e49f3151742b81945afe86b208ccd 100644 --- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml +++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml index cb815acd642f2981adb5943bb28968345d664e4f..2958bb531f203b02023eb2f38b6d36d3e46d996e 100644 --- a/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml +++ b/examples/EAGLE_ICs/EAGLE_25_low_res/eagle_25.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml index f102e15376882b75f329d8919ddbe57a433e834f..99120ad01d6c18ebfaeb074851b900e846cda6fc 100644 --- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml +++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml index a6bd67f534c792364af54ed5a5ba6b7393e96698..f6a3ec3049ddc9036e35fa28fecdbc22786f3bdb 100644 --- a/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml +++ b/examples/EAGLE_ICs/EAGLE_50_low_res/eagle_50.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml index a8f4ba314f94259bb2127f94509ad8ccc22ffd8b..5fdc200155d3912b7da12efc07ffbd1241225fb3 100644 --- a/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml +++ b/examples/EAGLE_ICs/EAGLE_6/eagle_6.yml @@ -30,6 +30,8 @@ Snapshots: output_list_on: 1 output_list: ./output_list.txt compression: 4 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_ICs/EAGLE_6/output_list.txt b/examples/EAGLE_ICs/EAGLE_6/output_list.txt index 151f609a9ad4dfbb142e28e721c1a706facdb511..592ab8483d015fe1bfafe5cc603fabc230b25589 100644 --- a/examples/EAGLE_ICs/EAGLE_6/output_list.txt +++ b/examples/EAGLE_ICs/EAGLE_6/output_list.txt @@ -1,5 +1,4 @@ # Redshift -126 18.08 15.28 13.06 diff --git a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml index 73adfed254f40866a37a65dacad614e4f7e27d52..14a42441f33462ab3ca698b35d69074010b11309 100644 --- a/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml +++ b/examples/EAGLE_low_z/EAGLE_100/eagle_100.yml @@ -34,6 +34,8 @@ Snapshots: scale_factor_first: 0.91 # Scale-factor of the first snaphot (cosmological run) time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units) delta_time: 1.01 # Time difference between consecutive outputs (in internal units) + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml index 156a2d3f459b68c1701aca364c5e3e2c7d95255f..24e6c7a53aae26663049cf838a930b626d6ab4f6 100644 --- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml +++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml @@ -32,6 +32,8 @@ Snapshots: time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units) delta_time: 1.01 # Time difference between consecutive outputs (in internal units) compression: 1 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot Scheduler: links_per_tasks: 500 diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml index fdf017e403a84e980fab7ea05c34b3ab9813b6f3..d79359fc684ce4e042c7827ccceb28c7e3d8db15 100644 --- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml +++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml @@ -42,6 +42,8 @@ Snapshots: scale_factor_first: 0.91 # Scale-factor of the first snaphot (cosmological run) time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units) delta_time: 1.01 # Time difference between consecutive outputs (in internal units) + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml index fc798ef4226343fd45f27f2a68681f21b040321c..c84621a62512f0429defe6bbbff18abd4296e281 100644 --- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml +++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml @@ -34,6 +34,8 @@ Snapshots: scale_factor_first: 0.91 # Scale-factor of the first snaphot (cosmological run) time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units) delta_time: 1.01 # Time difference between consecutive outputs (in internal units) + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml index db36cd7a0123a6222304e6dcda24a550bd035d1f..bc1b78866b0785edbb12871285f192989480c159 100644 --- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml +++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml @@ -45,6 +45,8 @@ Snapshots: delta_time: 1.01 # Time difference between consecutive outputs (in internal units) compression: 1 distributed: 1 + recording_triggers_part: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot + recording_triggers_bpart: [1.0227e-4, 1.0227e-5] # Recording starts 100M and 10M years before a snapshot # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml index 3c024595a00b84c5a896579b8724486d6571c346..c14d46ed5d5e434a6d18aa4023673a8d203c0077 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml @@ -30,9 +30,11 @@ TimeIntegration: Snapshots: basename: output # Common part of the name of output files time_first: 0. # Time of the first output if non-cosmological time-integration (in internal units) - delta_time: 0.001 # Time difference between consecutive outputs (in internal units) + delta_time: 0.001023 # Time difference between consecutive outputs (in internal units) 0.001023 TU = 1 Myr compression: 4 # Compress the snapshots - + recording_triggers_part: [-1, -1] # Not recording as we have many snapshots + recording_triggers_bpart: [-1, -1] # Not recording as we have many snapshots + # Parameters governing the conserved quantities statistics Statistics: delta_time: 1e-2 # Time between statistics output diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml index e1e522c3f749c522ae90f4a41e7c7941be9fa1f1..fb8936126da7993f5e46a2f6cf84efe5bbf0a313 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_starformation/isolated_galaxy.yml @@ -26,7 +26,10 @@ TimeIntegration: Snapshots: basename: output # Common part of the name of output files time_first: 0. # (Optional) Time of the first output if non-cosmological time-integration (in internal units) - delta_time: 0.001 # Time difference between consecutive outputs (in internal units) + delta_time: 0.001023 # Time difference between consecutive outputs (in internal units) 0.001023 TU = 1 Myr + compression: 4 # Compress the snapshots + recording_triggers_part: [-1, -1] # Not recording as we have many snapshots + recording_triggers_bpart: [-1, -1] # Not recording as we have many snapshots # Parameters governing the conserved quantities statistics Statistics: diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml index 0c71a23587eb3ba4cd707c621c1150b54b29d9c6..770353d87ef8242e59ca40cc80123cbfde64898f 100644 --- a/examples/parameter_example.yml +++ b/examples/parameter_example.yml @@ -188,7 +188,9 @@ Snapshots: select_output: selectoutput.yml # (Optional) File containing information to select outputs with (see documentation in the "Output Selection" section) run_on_dump: 0 # (Optional) Run the dump_command each time that a snapshot is dumped? dump_command: ./submit_velociraptor.sh # (Optional) Command to run each time that a snapshot is dumped. - + recording_triggers_part: [1e-3, 1e-2] # (Optional) Time before the snapshots where the trigger for gas particle tracers start (in internal units). + recording_triggers_spart: [1e-3, 1e-2] # (Optional) Time before the snapshots where the trigger for star particle tracers start (in internal units). + recording_triggers_bpart: [1e-3, 1e-2] # (Optional) Time before the snapshots where the trigger for BH particle tracers start (in internal units). # Parameters governing the CSDS snapshot system CSDS: diff --git a/src/Makefile.am b/src/Makefile.am index 43232da391459311f4c58dfd86a956f93e65482c..80a24fd178407b88df97c52b6fb4580a08a24abf 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -52,15 +52,17 @@ include_HEADERS += csds.h active.h timeline.h xmf.h gravity_properties.h gravity include_HEADERS += gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h include_HEADERS += chemistry.h chemistry_io.h chemistry_struct.h chemistry_debug.h cosmology.h restart.h space_getsid.h utilities.h include_HEADERS += cbrt.h exp10.h velociraptor_interface.h swift_velociraptor_part.h output_list.h -include_HEADERS += csds_io.h tracers_io.h tracers.h tracers_struct.h tracers_debug.h star_formation_io.h star_formation_debug.h extra_io.h +include_HEADERS += csds_io.h +include_HEADERS += tracers_io.h tracers.h tracers_triggers.h tracers_struct.h tracers_debug.h +include_HEADERS += star_formation_io.h star_formation_debug.h extra_io.h include_HEADERS += fof.h fof_struct.h fof_io.h fof_catalogue_io.h include_HEADERS += multipole.h multipole_accept.h multipole_struct.h binomial.h integer_power.h sincos.h include_HEADERS += star_formation_struct.h star_formation.h star_formation_iact.h include_HEADERS += star_formation_logger.h star_formation_logger_struct.h include_HEADERS += pressure_floor.h pressure_floor_struct.h pressure_floor_iact.h pressure_floor_debug.h include_HEADERS += velociraptor_struct.h velociraptor_io.h random.h memuse.h mpiuse.h memuse_rnodes.h -include_HEADERS += black_holes.h black_holes_io.h black_holes_properties.h black_holes_struct.h black_holes_debug.h -include_HEADERS += feedback.h feedback_struct.h feedback_properties.h feedback_debug.h +include_HEADERS += black_holes.h black_holes_iact.h black_holes_io.h black_holes_properties.h black_holes_struct.h black_holes_debug.h +include_HEADERS += feedback.h feedback_new_stars.h feedback_struct.h feedback_properties.h feedback_debug.h feedback_iact.h include_HEADERS += space_unique_id.h line_of_sight.h io_compression.h include_HEADERS += rays.h rays_struct.h include_HEADERS += sink.h sink_struct.h sink_io.h sink_properties.h sink_debug.h diff --git a/src/active.h b/src/active.h index bbdedf77d79bd9816849e4ce555a63c7e0833395..0556864feebabf0967657e256b0fbd55374b55ab 100644 --- a/src/active.h +++ b/src/active.h @@ -25,7 +25,7 @@ /* Local includes. */ #include "cell.h" #include "engine.h" -#include "feedback.h" +#include "feedback_new_stars.h" #include "part.h" #include "timeline.h" diff --git a/src/black_holes.h b/src/black_holes.h index 2c0f04ea7d3b7a7a846cb76fcc868c0c341fc87f..4051727192b99c512c79f5c0391eac4aa18857d9 100644 --- a/src/black_holes.h +++ b/src/black_holes.h @@ -25,15 +25,12 @@ /* Select the correct BH model */ #if defined(BLACK_HOLES_NONE) #include "./black_holes/Default/black_holes.h" -#include "./black_holes/Default/black_holes_iact.h" #elif defined(BLACK_HOLES_EAGLE) #include "./black_holes/EAGLE/black_holes.h" -#include "./black_holes/EAGLE/black_holes_iact.h" #elif defined(BLACK_HOLES_SPIN_JET) #include "./black_holes/SPIN_JET/black_holes.h" -#include "./black_holes/SPIN_JET/black_holes_iact.h" #else #error "Invalid choice of black hole model" #endif -#endif +#endif /* SWIFT_BLACK_HOLES_H */ diff --git a/src/black_holes/Default/black_holes_part.h b/src/black_holes/Default/black_holes_part.h index 777b13f8ccd78a8572e1b8ca066c85a30a6a830b..d8e379dcec878993c41c222ef7bc6de2c9e755ab 100644 --- a/src/black_holes/Default/black_holes_part.h +++ b/src/black_holes/Default/black_holes_part.h @@ -74,6 +74,9 @@ struct bpart { /*! Black holes merger information (e.g. merging ID) */ struct black_holes_bpart_data merger_data; + /*! Tracer structure */ + struct tracers_bpart_data tracers_data; + #ifdef SWIFT_DEBUG_CHECKS /* Time of the last drift */ diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h index cdff3735f6ffeea03f582882df650c635c8c9136..9a588be7eae825f4b9b7894755dfe08574defe31 100644 --- a/src/black_holes/EAGLE/black_holes_io.h +++ b/src/black_holes/EAGLE/black_holes_io.h @@ -155,7 +155,6 @@ INLINE static void convert_bpart_gas_temperatures(const struct engine* e, ret[0] = bp->internal_energy_gas * cosmo->a_factor_internal_energy / props->temp_to_u_factor; } - /** * @brief Specifies which b-particle fields to write to a dataset * @@ -167,7 +166,7 @@ INLINE static void convert_bpart_gas_temperatures(const struct engine* e, INLINE static void black_holes_write_particles(const struct bpart* bparts, struct io_props* list, int* num_fields, - int with_cosmology) { + const int with_cosmology) { /* Say how much we want to write */ *num_fields = 44; diff --git a/src/black_holes/EAGLE/black_holes_part.h b/src/black_holes/EAGLE/black_holes_part.h index ecb5541fdad20b8e3549ccfaf2b0f478849da09a..51bd44ec8b12efc4df8c9ebc614c8aa08a7eed78 100644 --- a/src/black_holes/EAGLE/black_holes_part.h +++ b/src/black_holes/EAGLE/black_holes_part.h @@ -274,6 +274,9 @@ struct bpart { /*! Isotropic AGN feedback information */ struct ray_data rays[eagle_blackhole_number_of_rays]; + /*! Tracer structure */ + struct tracers_bpart_data tracers_data; + #ifdef SWIFT_DEBUG_CHECKS /* Time of the last drift */ diff --git a/src/black_holes/SPIN_JET/black_holes_io.h b/src/black_holes/SPIN_JET/black_holes_io.h index 55d43c3305cfc0dfe2eee4afa256ca1dfe02e6c2..d3e80a34d0edb213c61d4794bb84208cb60939e2 100644 --- a/src/black_holes/SPIN_JET/black_holes_io.h +++ b/src/black_holes/SPIN_JET/black_holes_io.h @@ -173,7 +173,7 @@ INLINE static void convert_bpart_gas_velocity_curl(const struct engine* e, INLINE static void black_holes_write_particles(const struct bpart* bparts, struct io_props* list, int* num_fields, - int with_cosmology) { + const int with_cosmology) { /* Say how much we want to write */ *num_fields = 52; diff --git a/src/black_holes/SPIN_JET/black_holes_part.h b/src/black_holes/SPIN_JET/black_holes_part.h index c97a8ea105f753a2544d689729a8e879cc23eb1f..9df709e023e5350432ce59819ec8e612f2e454ff 100644 --- a/src/black_holes/SPIN_JET/black_holes_part.h +++ b/src/black_holes/SPIN_JET/black_holes_part.h @@ -331,6 +331,9 @@ struct bpart { oppposite direction from the spin axis. */ struct ray_data rays_jet_pos[spinjet_blackhole_number_of_rays]; + /*! Tracer structure */ + struct tracers_bpart_data tracers_data; + #ifdef SWIFT_DEBUG_CHECKS /* Time of the last drift */ diff --git a/src/black_holes/SPIN_JET/black_holes_spin.h b/src/black_holes/SPIN_JET/black_holes_spin.h index d1289686422944c165c9d49feaa79334d3d93fa3..29b5019c9fb873ce83423133aab1d7db8296a2df 100644 --- a/src/black_holes/SPIN_JET/black_holes_spin.h +++ b/src/black_holes/SPIN_JET/black_holes_spin.h @@ -860,9 +860,8 @@ __attribute__((always_inline)) INLINE static float black_hole_feedback_dv_jet( props->v_jet_BH_mass_scaling_slope); /* Get the critical density and virial overdensity at this redshift */ - const float critical_density = - cosmo->critical_density const float overdensity = - cosmo->overdensity_BN98; + const float critical_density = cosmo->critical_density; + const float overdensity = cosmo->overdensity_BN98; /* Gather the previous factors and compute the virial radius, virial velocity and finally the sound speed in the hot gas */ diff --git a/src/black_holes_iact.h b/src/black_holes_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..01fcd73160bb430167ce27c7745d82726e9be3d1 --- /dev/null +++ b/src/black_holes_iact.h @@ -0,0 +1,36 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2022 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_BLACK_HOLES_IACT_H +#define SWIFT_BLACK_HOLES_IACT_H + +/* Config parameters. */ +#include <config.h> + +/* Select the correct BH model */ +#if defined(BLACK_HOLES_NONE) +#include "./black_holes/Default/black_holes_iact.h" +#elif defined(BLACK_HOLES_EAGLE) +#include "./black_holes/EAGLE/black_holes_iact.h" +#elif defined(BLACK_HOLES_SPIN_JET) +#include "./black_holes/SPIN_JET/black_holes_iact.h" +#else +#error "Invalid choice of black hole model" +#endif + +#endif /* SWIFT_BLACK_HOLES_IACT_H */ diff --git a/src/common_io.c b/src/common_io.c index d07d5780251e412c05b98ada95701b651e2bd94e..e4db20223f0758c0deb4da88f6c6c6c157c791c3 100644 --- a/src/common_io.c +++ b/src/common_io.c @@ -1825,6 +1825,8 @@ void io_select_bh_fields(const struct bpart* const bparts, *num_fields += particle_splitting_write_bparticles(bparts, list + *num_fields); *num_fields += chemistry_write_bparticles(bparts, list + *num_fields); + *num_fields += + tracers_write_bparticles(bparts, list + *num_fields, with_cosmology); if (with_fof) { *num_fields += fof_write_bparts(bparts, list + *num_fields); } diff --git a/src/engine.c b/src/engine.c index 68e10d3805abf4808ed89b47a8f131fe87785455..994a29e4aed181c5b4992136804c0cd74a131268 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2347,6 +2347,9 @@ int engine_step(struct engine *e) { hydro_props_update(e->hydro_properties, e->gravity_properties, e->cosmology); + /* Check for any snapshot triggers */ + engine_io_check_snapshot_triggers(e); + if (e->verbose) message("Updating general quantities took %.3f %s", clocks_from_ticks(getticks() - tic_updates), clocks_getunit()); @@ -3066,6 +3069,18 @@ void engine_init( e->min_active_bin_subcycle = 1; e->internal_units = internal_units; e->output_list_snapshots = NULL; + if (num_snapshot_triggers_part) + parser_get_param_double_array(params, "Snapshots:recording_triggers_part", + num_snapshot_triggers_part, + e->snapshot_recording_triggers_part); + if (num_snapshot_triggers_spart) + parser_get_param_double_array(params, "Snapshots:recording_triggers_spart", + num_snapshot_triggers_spart, + e->snapshot_recording_triggers_spart); + if (num_snapshot_triggers_bpart) + parser_get_param_double_array(params, "Snapshots:recording_triggers_bpart", + num_snapshot_triggers_bpart, + e->snapshot_recording_triggers_bpart); e->a_first_snapshot = parser_get_opt_param_double(params, "Snapshots:scale_factor_first", 0.1); e->time_first_snapshot = diff --git a/src/engine.h b/src/engine.h index 5fa8db8508b351276b446ee8ff0055636bd93f52..d6ab7095cc091ddc5f399f6630e57fc3e4f8c716 100644 --- a/src/engine.h +++ b/src/engine.h @@ -47,6 +47,7 @@ #include "scheduler.h" #include "space.h" #include "task.h" +#include "tracers_triggers.h" #include "units.h" #include "velociraptor_interface.h" @@ -355,6 +356,16 @@ struct engine { double snapshot_delta_from_edge; int snapshot_output_count; + /* Snapshot recording trigger mechanism counters */ + double snapshot_recording_triggers_part[num_snapshot_triggers_part]; + int snapshot_recording_triggers_started_part[num_snapshot_triggers_part]; + + double snapshot_recording_triggers_spart[num_snapshot_triggers_spart]; + int snapshot_recording_triggers_started_spart[num_snapshot_triggers_spart]; + + double snapshot_recording_triggers_bpart[num_snapshot_triggers_bpart]; + int snapshot_recording_triggers_started_bpart[num_snapshot_triggers_bpart]; + /* Metadata from the ICs */ struct ic_info *ics_metadata; @@ -675,6 +686,8 @@ void engine_reconstruct_multipoles(struct engine *e); void engine_allocate_foreign_particles(struct engine *e, const int fof); void engine_print_stats(struct engine *e); void engine_io(struct engine *e); +void engine_io_check_snapshot_triggers(struct engine *e); +void engine_io_verify_triggers_size(const struct engine *e); void engine_collect_end_of_step(struct engine *e, int apply); void engine_collect_end_of_sub_cycle(struct engine *e); void engine_dump_snapshot(struct engine *e); diff --git a/src/engine_config.c b/src/engine_config.c index e6b028fdee2ab2bbbea1a29b77d268a4def3095c..f47594101c60e3efddb1d6969bf705bd81a6c650 100644 --- a/src/engine_config.c +++ b/src/engine_config.c @@ -613,6 +613,10 @@ void engine_config(int restart, int fof, struct engine *e, /* Read (or re-read the list of outputs */ engine_init_output_lists(e, params, e->output_options); + /* Verify that the snapshot triggers are compatible with the + * output list */ + engine_io_verify_triggers_size(e); + /* Check whether output quantities make sense */ if (e->policy & engine_policy_cosmology) { diff --git a/src/engine_io.c b/src/engine_io.c index 5b934ce3502e4b0fa9cd1d202c464b88282fd138..b45e8fd24886fa4fc72f7322341b1aad5d54c3ae 100644 --- a/src/engine_io.c +++ b/src/engine_io.c @@ -31,6 +31,7 @@ #include "engine.h" /* Local headers. */ +#include "active.h" #include "csds_io.h" #include "distributed_io.h" #include "kick.h" @@ -41,9 +42,128 @@ #include "power_spectrum.h" #include "serial_io.h" #include "single_io.h" +#include "tracers.h" +/* Standard includes */ #include <stdio.h> +/** + * @brief Finalize the quantities recorded via the snapshot triggers + * + * This adds the small amount of time between the particles' last start + * of step and the current (snapshot) time. + * + * @param e The #engine to act on. + */ +void engine_finalize_trigger_recordings(struct engine *e) { + + const int with_cosmology = (e->policy & engine_policy_cosmology); + struct space *s = e->s; + + /* Finish the recording period for part triggers */ + if (num_snapshot_triggers_part) { + for (size_t k = 0; k < s->nr_parts; ++k) { + + /* Get a handle on the part. */ + struct part *p = &s->parts[k]; + struct xpart *xp = &s->xparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, p->time_bin); + + /* Escape inhibited particles */ + if (part_is_inhibited(p, e)) continue; + + /* We need to escape the special case of a particle that + * actually ended its time-step on this very step */ + if (e->ti_current - ti_begin == get_integer_timestep(p->time_bin)) + continue; + + /* Time from the start of the particle's step to the snapshot (aka. + * current time) */ + double missing_time; + if (with_cosmology) { + missing_time = + cosmology_get_delta_time(e->cosmology, ti_begin, e->ti_current); + } else { + missing_time = (e->ti_current - ti_begin) * e->time_base; + } + + tracers_after_timestep_part( + p, xp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, e->hydro_properties, e->cooling_func, e->time, + missing_time, e->snapshot_recording_triggers_started_part); + } + } + + /* Finish the recording period for spart triggers */ + if (num_snapshot_triggers_spart) { + for (size_t k = 0; k < s->nr_sparts; ++k) { + + /* Get a handle on the part. */ + struct spart *sp = &s->sparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, sp->time_bin); + + /* Escape inhibited particles */ + if (spart_is_inhibited(sp, e)) continue; + + /* We need to escape the special case of a particle that + * actually ended its time-step on this very step */ + if (e->ti_current - ti_begin == get_integer_timestep(sp->time_bin)) + continue; + + /* Time from the start of the particle's step to the snapshot (aka. + * current time) */ + double missing_time; + if (with_cosmology) { + missing_time = + cosmology_get_delta_time(e->cosmology, ti_begin, e->ti_current); + } else { + missing_time = (e->ti_current - ti_begin) * e->time_base; + } + + tracers_after_timestep_spart( + sp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, missing_time, + e->snapshot_recording_triggers_started_spart); + } + } + + /* Finish the recording period for spart triggers */ + if (num_snapshot_triggers_bpart) { + for (size_t k = 0; k < s->nr_bparts; ++k) { + + /* Get a handle on the part. */ + struct bpart *bp = &s->bparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, bp->time_bin); + + /* Escape inhibited particles */ + if (bpart_is_inhibited(bp, e)) continue; + + /* We need to escape the special case of a particle that + * actually ended its time-step on this very step */ + if (e->ti_current - ti_begin == get_integer_timestep(bp->time_bin)) + continue; + + /* Time from the start of the particle's step to the snapshot (aka. + * current time) */ + double missing_time; + if (with_cosmology) { + missing_time = + cosmology_get_delta_time(e->cosmology, ti_begin, e->ti_current); + } else { + missing_time = (e->ti_current - ti_begin) * e->time_base; + } + + tracers_after_timestep_bpart( + bp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, missing_time, + e->snapshot_recording_triggers_started_bpart); + } + } +} + /** * @brief dump restart files if it is time to do so and dumps are enabled. * @@ -184,9 +304,12 @@ void engine_dump_snapshot(struct engine *e) { engine_collect_stars_counter(e); #endif + /* Finalize the recording periods */ + engine_finalize_trigger_recordings(e); + /* Get time-step since the last mesh kick */ if ((e->policy & engine_policy_self_gravity) && e->s->periodic) { - const int with_cosmology = e->policy & engine_policy_cosmology; + const int with_cosmology = (e->policy & engine_policy_cosmology); e->dt_kick_grav_mesh_for_io = kick_get_grav_kick_dt(e->mesh->ti_beg_mesh_next, e->ti_current, @@ -220,6 +343,22 @@ void engine_dump_snapshot(struct engine *e) { #endif #endif + /* Cancel any triggers that are switched on */ + if (num_snapshot_triggers_part || num_snapshot_triggers_spart || + num_snapshot_triggers_bpart) { + + /* Reset the trigger flags */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) + e->snapshot_recording_triggers_started_part[i] = 0; + for (int i = 0; i < num_snapshot_triggers_spart; ++i) + e->snapshot_recording_triggers_started_spart[i] = 0; + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) + e->snapshot_recording_triggers_started_bpart[i] = 0; + + /* Reser the tracers themselves */ + space_after_snap_tracer(e->s, e->verbose); + } + /* Flag that we dumped a snapshot */ e->step_props |= engine_step_prop_snapshot; @@ -1059,3 +1198,352 @@ void engine_init_output_lists(struct engine *e, struct swift_params *params, } } } + +/** + * @brief Checks whether we passed a certain delta time before the next snapshot + * and need to trigger a recording. + * + * If a recording has to start, we also loop over the particles to correct for + * the time between the particles' end of time-step and the actual start of + * trigger. + * + * @param e The #engine. + */ +void engine_io_check_snapshot_triggers(struct engine *e) { + + struct space *s = e->s; + const int with_cosmology = (e->policy & engine_policy_cosmology); + + /* Time until the next snapshot */ + double time_to_next_snap; + if (e->policy & engine_policy_cosmology) { + time_to_next_snap = cosmology_get_delta_time(e->cosmology, e->ti_current, + e->ti_next_snapshot); + } else { + time_to_next_snap = (e->ti_next_snapshot - e->ti_current) * e->time_base; + } + + /* Should any not yet switched on trigger be activated? (part version) */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + + if (time_to_next_snap <= e->snapshot_recording_triggers_part[i] && + e->snapshot_recording_triggers_part[i] > 0. && + !e->snapshot_recording_triggers_started_part[i]) { + e->snapshot_recording_triggers_started_part[i] = 1; + + /* Be vocal about this */ + if (e->verbose) + message( + "Snapshot will be dumped in %e U_t. Recording trigger for part " + "activated.", + e->snapshot_recording_triggers_part[i]); + + /* We now need to loop over the particles to preemptively deduct the extra + * time logged between the particles' start of step and the actual start + * of the trigger */ + for (size_t k = 0; k < s->nr_parts; ++k) { + + /* Get a handle on the part. */ + struct part *p = &s->parts[k]; + struct xpart *xp = &s->xparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, p->time_bin); + + /* Escape inhibited particles */ + if (part_is_inhibited(p, e)) continue; + + /* Time from the start of the particle's step to the snapshot */ + double total_time; + if (with_cosmology) { + total_time = cosmology_get_delta_time(e->cosmology, ti_begin, + e->ti_next_snapshot); + } else { + total_time = (e->ti_next_snapshot - ti_begin) * e->time_base; + } + + /* Time to deduct = time since the start of the step - trigger time */ + const double time_to_remove = + total_time - e->snapshot_recording_triggers_part[i]; + +#ifdef SWIFT_DEBUG_CHECKS + if (time_to_remove < 0.) + error("Invalid time to deduct! %e", time_to_remove); +#endif + + /* Note that we need to use a separate array (not the raw + * e->snapshot_recording_triggers_part) as we only want to + * update one entry */ + int my_temp_array[num_snapshot_triggers_part]; + memset(my_temp_array, 0, sizeof(int) * num_snapshot_triggers_part); + my_temp_array[i] = 1; + + tracers_after_timestep_part( + p, xp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, e->hydro_properties, e->cooling_func, e->time, + -time_to_remove, my_temp_array); + } + } + } + + /* Should any not yet switched on trigger be activated? (spart version) */ + for (int i = 0; i < num_snapshot_triggers_spart; ++i) { + + if (time_to_next_snap <= e->snapshot_recording_triggers_spart[i] && + e->snapshot_recording_triggers_spart[i] > 0. && + !e->snapshot_recording_triggers_started_spart[i]) { + e->snapshot_recording_triggers_started_spart[i] = 1; + + /* Be vocal about this */ + if (e->verbose) + message( + "Snapshot will be dumped in %e U_t. Recording trigger for spart " + "activated.", + e->snapshot_recording_triggers_spart[i]); + + /* We now need to loop over the particles to preemptively deduct the extra + * time logged between the particles' start of step and the actual start + * of the trigger */ + for (size_t k = 0; k < s->nr_sparts; ++k) { + + /* Get a handle on the spart. */ + struct spart *sp = &s->sparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, sp->time_bin); + + /* Escape inhibited particles */ + if (spart_is_inhibited(sp, e)) continue; + + /* Time from the start of the particle's step to the snapshot */ + double total_time; + if (with_cosmology) { + total_time = cosmology_get_delta_time(e->cosmology, ti_begin, + e->ti_next_snapshot); + } else { + total_time = (e->ti_next_snapshot - ti_begin) * e->time_base; + } + + /* Time to deduct = time since the start of the step - trigger time */ + const double time_to_remove = + total_time - e->snapshot_recording_triggers_part[i]; + +#ifdef SWIFT_DEBUG_CHECKS + if (time_to_remove < 0.) + error("Invalid time to deduct! %e", time_to_remove); +#endif + + /* Note that we need to use a separate array (not the raw + * e->snapshot_recording_triggers_part) as we only want to + * update one entry */ + int my_temp_array[num_snapshot_triggers_spart]; + memset(my_temp_array, 0, sizeof(int) * num_snapshot_triggers_spart); + my_temp_array[i] = 1; + + tracers_after_timestep_spart( + sp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, -time_to_remove, my_temp_array); + } + } + } + + /* Should any not yet switched on trigger be activated? (bpart version) */ + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) { + + if (time_to_next_snap <= e->snapshot_recording_triggers_bpart[i] && + e->snapshot_recording_triggers_bpart[i] > 0. && + !e->snapshot_recording_triggers_started_bpart[i]) { + e->snapshot_recording_triggers_started_bpart[i] = 1; + + /* Be vocal about this */ + if (e->verbose) + message( + "Snapshot will be dumped in %e U_t. Recording trigger for bpart " + "activated.", + e->snapshot_recording_triggers_bpart[i]); + + /* We now need to loop over the particles to preemptively deduct the extra + * time logged between the particles' start of step and the actual start + * of the trigger */ + for (size_t k = 0; k < s->nr_bparts; ++k) { + + /* Get a handle on the bpart. */ + struct bpart *bp = &s->bparts[k]; + const integertime_t ti_begin = + get_integer_time_begin(e->ti_current, bp->time_bin); + + /* Escape inhibited particles */ + if (bpart_is_inhibited(bp, e)) continue; + + /* Time from the start of the particle's step to the snapshot */ + double total_time; + if (with_cosmology) { + total_time = cosmology_get_delta_time(e->cosmology, ti_begin, + e->ti_next_snapshot); + } else { + total_time = (e->ti_next_snapshot - ti_begin) * e->time_base; + } + + /* Time to deduct = time since the start of the step - trigger time */ + const double time_to_remove = + total_time - e->snapshot_recording_triggers_part[i]; + +#ifdef SWIFT_DEBUG_CHECKS + if (time_to_remove < 0.) + error("Invalid time to deduct! %e", time_to_remove); +#endif + + /* Note that we need to use a separate array (not the raw + * e->snapshot_recording_triggers_part) as we only want to + * update one entry */ + int my_temp_array[num_snapshot_triggers_bpart]; + memset(my_temp_array, 0, sizeof(int) * num_snapshot_triggers_bpart); + my_temp_array[i] = 1; + + tracers_after_timestep_bpart( + bp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, -time_to_remove, my_temp_array); + } + } + } +} + +void check_triggers(const struct engine *e, const double snapshot_delta) { + + /* Verify the part triggers */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + if (e->snapshot_recording_triggers_part[i] > snapshot_delta) { + error( + "The recording time for the part trigger %d (dt=%e) is longer " + "than the time between snapshots (dt=%e)", + i, e->snapshot_recording_triggers_part[i], snapshot_delta); + } + } + + /* Verify the spart triggers */ + for (int i = 0; i < num_snapshot_triggers_spart; ++i) { + if (e->snapshot_recording_triggers_spart[i] > snapshot_delta) { + error( + "The recording time for the spart trigger %d (dt=%e) is longer " + "than the time between snapshots (dt=%e)", + i, e->snapshot_recording_triggers_spart[i], snapshot_delta); + } + } + + /* Verify the bpart triggers */ + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) { + if (e->snapshot_recording_triggers_bpart[i] > snapshot_delta) { + error( + "The recording time for the bpart trigger %d (dt=%e) is longer " + "than the time between snapshots (dt=%e)", + i, e->snapshot_recording_triggers_bpart[i], snapshot_delta); + } + } +} + +/** + * @brief Verify that the triggers provided by the user are not + * overlapping with the snapshot times. + * + * Deals witht the case of a regular snapshot pattern and output list. + * + * @param e The #engine. + */ +void engine_io_verify_triggers_size(const struct engine *e) { + + /* Stop if we have nothing to do */ + if (num_snapshot_triggers_part == 0 && num_snapshot_triggers_spart == 0 && + num_snapshot_triggers_bpart == 0) + return; + + if (e->output_list_snapshots) { + + /* Case where an output list is used */ + + if (e->policy & engine_policy_cosmology) { + + /* Loop over all the snapshot gaps */ + for (size_t i = 0; i < e->output_list_snapshots->size - 1; ++i) { + + const double a = e->output_list_snapshots->times[i]; + const double a_next = e->output_list_snapshots->times[i + 1]; + + const double snapshot_delta = + cosmology_get_delta_time_from_scale_factors(e->cosmology, a, + a_next); + /* Do the verification */ + check_triggers(e, snapshot_delta); + } + + /* Also verify what happens before the first snapshot */ + const double snapshot_delta_first = + cosmology_get_delta_time_from_scale_factors( + e->cosmology, e->time_begin, e->output_list_snapshots->times[0]); + + /* Do the verification */ + check_triggers(e, snapshot_delta_first); + + } else { + + /* Loop over all the snapshot gaps */ + for (size_t i = 0; i < e->output_list_snapshots->size - 1; ++i) { + + const double time = e->output_list_snapshots->times[i]; + const double time_next = e->output_list_snapshots->times[i + 1]; + + const double snapshot_delta = time_next - time; + + /* Do the verification */ + check_triggers(e, snapshot_delta); + } + + /* Also verify what happens before the first snapshot */ + const double snapshot_delta_first = + e->output_list_snapshots->times[0] - e->time_begin; + + /* Do the verification */ + check_triggers(e, snapshot_delta_first); + } + + } else { + + /* Case where we have a regular pattern */ + + if (e->policy & engine_policy_cosmology) { + + /* Loop over all the snapshot gaps */ + for (double a = e->a_first_snapshot; a <= e->time_end; + a *= e->delta_time_snapshot) { + + const double a_next = a * e->delta_time_snapshot; + const double snapshot_delta = + cosmology_get_delta_time_from_scale_factors(e->cosmology, a, + a_next); + + /* Do the verification */ + check_triggers(e, snapshot_delta); + } + + /* Also verify what happens before the first snapshot */ + const double snapshot_delta_first = + cosmology_get_delta_time_from_scale_factors( + e->cosmology, e->time_begin, e->a_first_snapshot); + + /* Do the verification */ + check_triggers(e, snapshot_delta_first); + + } else { + + /* Simplest case: the delta is constant */ + const double snapshot_delta = e->delta_time_snapshot; + + /* Do the verification */ + check_triggers(e, snapshot_delta); + + /* Also verify what happens before the first snapshot */ + const double snapshot_delta_first = + e->time_first_snapshot - e->time_begin; + + /* Do the verification */ + check_triggers(e, snapshot_delta_first); + } + } +} diff --git a/src/feedback.h b/src/feedback.h index 2cf74b61c7a21abf3355edf8818ea4441313c184..a403cd25a062d92ec9e5704b52ee39363f8794d0 100644 --- a/src/feedback.h +++ b/src/feedback.h @@ -25,21 +25,13 @@ /* Select the correct feedback model */ #if defined(FEEDBACK_NONE) #include "./feedback/none/feedback.h" -#include "./feedback/none/feedback_iact.h" -#define feedback_use_newborn_stars 0 #elif defined(FEEDBACK_EAGLE_THERMAL) #include "./feedback/EAGLE_thermal/feedback.h" -#include "./feedback/EAGLE_thermal/feedback_iact.h" -#define feedback_use_newborn_stars 0 #elif defined(FEEDBACK_EAGLE_KINETIC) #include "./feedback/EAGLE_kinetic/feedback.h" -#include "./feedback/EAGLE_kinetic/feedback_iact.h" -#define feedback_use_newborn_stars 0 #define EXTRA_STAR_LOOPS #elif defined(FEEDBACK_GEAR) #include "./feedback/GEAR/feedback.h" -#include "./feedback/GEAR/feedback_iact.h" -#define feedback_use_newborn_stars 1 #else #error "Invalid choice of feedback model" #endif diff --git a/src/feedback/EAGLE_kinetic/feedback_iact.h b/src/feedback/EAGLE_kinetic/feedback_iact.h index b470572f584203fc28148ab30c07fd1dee89f39d..9efe800fbf421637f1f8b302f4b15ad3352270ed 100644 --- a/src/feedback/EAGLE_kinetic/feedback_iact.h +++ b/src/feedback/EAGLE_kinetic/feedback_iact.h @@ -20,6 +20,7 @@ #define SWIFT_EAGLE_FEEDBACK_KINETIC_IACT_H /* Local includes */ +#include "feedback.h" #include "random.h" #include "rays.h" #include "timestep_sync_part.h" diff --git a/src/feedback/EAGLE_thermal/feedback_iact.h b/src/feedback/EAGLE_thermal/feedback_iact.h index d1b84ca992575a79feae3970400dbc092d6bf987..89c91a63984be432776494b35e0794522da1246c 100644 --- a/src/feedback/EAGLE_thermal/feedback_iact.h +++ b/src/feedback/EAGLE_thermal/feedback_iact.h @@ -20,6 +20,7 @@ #define SWIFT_EAGLE_FEEDBACK_IACT_THERMAL_H /* Local includes */ +#include "feedback.h" #include "random.h" #include "rays.h" #include "timestep_sync_part.h" diff --git a/src/feedback/GEAR/feedback_iact.h b/src/feedback/GEAR/feedback_iact.h index 36f42a044e4b5941b233966c423e53babbb03f86..80c205d3f56c58efdf9a2e9d81400bf85e794c21 100644 --- a/src/feedback/GEAR/feedback_iact.h +++ b/src/feedback/GEAR/feedback_iact.h @@ -20,6 +20,7 @@ #define SWIFT_GEAR_FEEDBACK_IACT_H /* Local includes */ +#include "feedback.h" #include "hydro.h" #include "random.h" #include "timestep_sync_part.h" diff --git a/src/feedback/none/feedback_iact.h b/src/feedback/none/feedback_iact.h index 02041ad67ec1891b5c461b8313a20a4db064c37f..e630ff68048040d78f010365cb6e68761c060b4c 100644 --- a/src/feedback/none/feedback_iact.h +++ b/src/feedback/none/feedback_iact.h @@ -19,6 +19,9 @@ #ifndef SWIFT_NONE_FEEDBACK_IACT_H #define SWIFT_NONE_FEEDBACK_IACT_H +/* Local includes */ +#include "feedback.h" + /** * @brief Density interaction between two particles (non-symmetric). * diff --git a/src/feedback_iact.h b/src/feedback_iact.h new file mode 100644 index 0000000000000000000000000000000000000000..a6144ab2c58402d4feaeb1b3a8e53bd652c8be3f --- /dev/null +++ b/src/feedback_iact.h @@ -0,0 +1,38 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2022 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_FEEDBACK_IACT_H +#define SWIFT_FEEDBACK_IACT_H + +/* Config parameters. */ +#include <config.h> + +/* Select the correct feedback model */ +#if defined(FEEDBACK_NONE) +#include "./feedback/none/feedback_iact.h" +#elif defined(FEEDBACK_EAGLE_THERMAL) +#include "./feedback/EAGLE_thermal/feedback_iact.h" +#elif defined(FEEDBACK_EAGLE_KINETIC) +#include "./feedback/EAGLE_kinetic/feedback_iact.h" +#elif defined(FEEDBACK_GEAR) +#include "./feedback/GEAR/feedback_iact.h" +#else +#error "Invalid choice of feedback model" +#endif + +#endif /* SWIFT_FEEDBACK_IACT_H */ diff --git a/src/feedback_new_stars.h b/src/feedback_new_stars.h new file mode 100644 index 0000000000000000000000000000000000000000..010a32f2f74152e739b4299bcb96f05a4aa84d4b --- /dev/null +++ b/src/feedback_new_stars.h @@ -0,0 +1,39 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2022 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_FEEDBACK_NEW_STARS_H +#define SWIFT_FEEDBACK_NEW_STARS_H + +/* Config parameters. */ +#include <config.h> + +/* Select the correct feedback model and switches on/off the + * feedback from stars born in that very step. */ +#if defined(FEEDBACK_NONE) +#define feedback_use_newborn_stars 0 +#elif defined(FEEDBACK_EAGLE_THERMAL) +#define feedback_use_newborn_stars 0 +#elif defined(FEEDBACK_EAGLE_KINETIC) +#define feedback_use_newborn_stars 0 +#elif defined(FEEDBACK_GEAR) +#define feedback_use_newborn_stars 1 +#else +#error "Invalid choice of feedback model" +#endif + +#endif diff --git a/src/fof.c b/src/fof.c index de69efda987a01a804839216186ce0e4c9f8230d..c5da7e068af4ca1528316687d07ea43aa63ced8b 100644 --- a/src/fof.c +++ b/src/fof.c @@ -45,6 +45,7 @@ #include "proxy.h" #include "threadpool.h" #include "tools.h" +#include "tracers.h" #define fof_props_default_group_id 2147483647 #define fof_props_default_group_id_offset 1 @@ -2215,6 +2216,8 @@ void fof_seed_black_holes(const struct fof_props *props, /* Copy over all the gas properties that we want */ black_holes_create_from_gas(bp, bh_props, constants, cosmo, p, xp, s->e->ti_current); + tracers_first_init_bpart(bp, s->e->internal_units, + s->e->physical_constants, cosmo); /* Move to the next BH slot */ k++; diff --git a/src/runner_doiact_black_holes.c b/src/runner_doiact_black_holes.c index bed8b293d78dafb5608e15ae0e09ce07725fed1c..ac0414302407cd1c8cd436687dd408ce6432bce4 100644 --- a/src/runner_doiact_black_holes.c +++ b/src/runner_doiact_black_holes.c @@ -24,7 +24,7 @@ /* Local headers. */ #include "active.h" -#include "black_holes.h" +#include "black_holes_iact.h" #include "cell.h" #include "engine.h" #include "runner.h" diff --git a/src/runner_doiact_stars.c b/src/runner_doiact_stars.c index 00dd94cc24c79a4319810a7a619d26874cbc72f8..4720a96fbd02e9d93eb4ddb8b59c82157a3fdaf0 100644 --- a/src/runner_doiact_stars.c +++ b/src/runner_doiact_stars.c @@ -26,7 +26,7 @@ #include "active.h" #include "cell.h" #include "engine.h" -#include "feedback.h" +#include "feedback_iact.h" #include "rt.h" #include "runner.h" #include "space_getsid.h" diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c index 34c0eccb639e294df1ec9e37792434d32138f28b..47164fa64764a7e565a28cef9abbee364dc2537a 100644 --- a/src/runner_time_integration.c +++ b/src/runner_time_integration.c @@ -716,6 +716,15 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { error("Computing RT time-step of rogue particle"); } #endif + /* Old time-step length in physical units */ + const integertime_t ti_old_step = get_integer_timestep(p->time_bin); + double old_time_step_length; + if (with_cosmology) { + old_time_step_length = cosmology_get_delta_time( + e->cosmology, e->ti_current - ti_old_step, e->ti_current); + } else { + old_time_step_length = get_timestep(p->time_bin, e->time_base); + } /* Get new time-step */ const integertime_t ti_new_step = get_part_timestep(p, xp, e); @@ -726,13 +735,13 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { /* Update particle */ p->time_bin = get_time_bin(ti_new_step); - if (p->gpart != NULL) p->gpart->time_bin = p->time_bin; /* Update the tracers properties */ - tracers_after_timestep(p, xp, e->internal_units, e->physical_constants, - with_cosmology, e->cosmology, - e->hydro_properties, e->cooling_func, e->time); + tracers_after_timestep_part( + p, xp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, e->hydro_properties, e->cooling_func, e->time, + old_time_step_length, e->snapshot_recording_triggers_started_part); /* Number of updated particles */ updated++; @@ -903,6 +912,16 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { if (ti_end != ti_current) error("Computing time-step of rogue particle."); #endif + /* Old time-step length in physical units */ + const integertime_t ti_old_step = get_integer_timestep(sp->time_bin); + double old_time_step_length; + if (with_cosmology) { + old_time_step_length = cosmology_get_delta_time( + e->cosmology, e->ti_current - ti_old_step, e->ti_current); + } else { + old_time_step_length = get_timestep(sp->time_bin, e->time_base); + } + /* Get new time-step */ const integertime_t ti_new_step = get_spart_timestep(sp, e); @@ -910,6 +929,12 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { sp->time_bin = get_time_bin(ti_new_step); sp->gpart->time_bin = get_time_bin(ti_new_step); + /* Update the tracers properties */ + tracers_after_timestep_spart( + sp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, old_time_step_length, + e->snapshot_recording_triggers_started_spart); + /* Update feedback related counters */ if (with_feedback) { @@ -1031,6 +1056,16 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { if (ti_end != ti_current) error("Computing time-step of rogue particle."); #endif + /* Old time-step length in physical units */ + const integertime_t ti_old_step = get_integer_timestep(bp->time_bin); + double old_time_step_length; + if (with_cosmology) { + old_time_step_length = cosmology_get_delta_time( + e->cosmology, e->ti_current - ti_old_step, e->ti_current); + } else { + old_time_step_length = get_timestep(bp->time_bin, e->time_base); + } + /* Get new time-step */ const integertime_t ti_new_step = get_bpart_timestep(bp, e); @@ -1038,6 +1073,12 @@ void runner_do_timestep(struct runner *r, struct cell *c, const int timer) { bp->time_bin = get_time_bin(ti_new_step); bp->gpart->time_bin = get_time_bin(ti_new_step); + /* Update the tracers properties */ + tracers_after_timestep_bpart( + bp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, old_time_step_length, + e->snapshot_recording_triggers_started_bpart); + /* Number of updated s-particles */ b_updated++; g_updated++; @@ -1528,14 +1569,25 @@ void runner_do_sync(struct runner *r, struct cell *c, int force, new_time_bin = min(new_time_bin, e->max_active_bin); ti_new_step = get_integer_timestep(new_time_bin); + /* Time-step length in physical units */ + // MATTHIEU: TODO: think about this one! + double time_step_length; + if (with_cosmology) { + time_step_length = cosmology_get_delta_time( + e->cosmology, e->ti_current, e->ti_current + ti_new_step); + } else { + time_step_length = get_timestep(new_time_bin, e->time_base); + } + /* Update particle */ p->time_bin = new_time_bin; if (p->gpart != NULL) p->gpart->time_bin = new_time_bin; /* Update the tracers properties */ - tracers_after_timestep(p, xp, e->internal_units, e->physical_constants, - with_cosmology, e->cosmology, - e->hydro_properties, e->cooling_func, e->time); + tracers_after_timestep_part( + p, xp, e->internal_units, e->physical_constants, with_cosmology, + e->cosmology, e->hydro_properties, e->cooling_func, e->time, + 0 * time_step_length, e->snapshot_recording_triggers_started_part); #ifdef SWIFT_HYDRO_DENSITY_CHECKS p->limited_part = 1; diff --git a/src/space.c b/src/space.c index eadac16762523a983c608ffcc07e210ba0cc6abe..d7cf8aed7963e9374c4997e6466c30a474fb5dd1 100644 --- a/src/space.c +++ b/src/space.c @@ -42,6 +42,7 @@ /* Local headers. */ #include "active.h" #include "atomic.h" +#include "black_holes.h" #include "const.h" #include "cooling.h" #include "engine.h" @@ -59,6 +60,7 @@ #include "stars.h" #include "threadpool.h" #include "tools.h" +#include "tracers.h" /* Split size. */ int space_splitsize = space_splitsize_default; @@ -2319,6 +2321,23 @@ void space_reset_task_counters(struct space *s) { #endif } +/** + * @brief Call the post-snapshot tracer on all the particles. + * + * @param s The #space. + */ +void space_after_snap_tracer(struct space *s, int verbose) { + for (size_t i = 0; i < s->nr_parts; ++i) { + tracers_after_snapshot_part(&s->parts[i], &s->xparts[i]); + } + for (size_t i = 0; i < s->nr_sparts; ++i) { + tracers_after_snapshot_spart(&s->sparts[i]); + } + for (size_t i = 0; i < s->nr_bparts; ++i) { + tracers_after_snapshot_bpart(&s->bparts[i]); + } +} + /** * @brief Frees up the memory allocated for this #space */ diff --git a/src/space.h b/src/space.h index e4f6ad8a9164ad77aa65e5090ebd48ca213177c7..371c3427d5903fb33318e28bf7a2250751329351 100644 --- a/src/space.h +++ b/src/space.h @@ -422,6 +422,7 @@ void space_init_gparts(struct space *s, int verbose); void space_init_sparts(struct space *s, int verbose); void space_init_bparts(struct space *s, int verbose); void space_init_sinks(struct space *s, int verbose); +void space_after_snap_tracer(struct space *s, int verbose); void space_convert_quantities(struct space *s, int verbose); void space_convert_rt_quantities(struct space *s, int verbose); void space_link_cleanup(struct space *s); diff --git a/src/space_first_init.c b/src/space_first_init.c index 5e76749e54d5b537c5127d21f942cf916130fe8b..71dbd26b50d9e9c2e6a19603e6d19ea4f22d82b5 100644 --- a/src/space_first_init.c +++ b/src/space_first_init.c @@ -269,6 +269,8 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count, const struct cosmology *cosmo = e->cosmology; const struct stars_props *stars_properties = e->stars_properties; const struct feedback_props *feedback_properties = e->feedback_props; + const struct phys_const *phys_const = s->e->physical_constants; + const struct unit_system *us = s->e->internal_units; const float a_factor_vel = cosmo->a; /* Convert velocities to internal units */ @@ -311,6 +313,9 @@ void space_first_init_sparts_mapper(void *restrict map_data, int count, csds_part_data_init(&sp[k].csds_data); #endif + /* And the tracers */ + tracers_first_init_spart(&sp[k], us, phys_const, cosmo); + /* Also initialise the chemistry */ chemistry_first_init_spart(chemistry, &sp[k]); @@ -366,6 +371,8 @@ void space_first_init_bparts_mapper(void *restrict map_data, int count, const float initial_h = s->initial_bpart_h; const struct cosmology *cosmo = e->cosmology; + const struct phys_const *phys_const = s->e->physical_constants; + const struct unit_system *us = s->e->internal_units; const float a_factor_vel = cosmo->a; /* Convert velocities to internal units */ @@ -403,6 +410,9 @@ void space_first_init_bparts_mapper(void *restrict map_data, int count, black_holes_first_init_bpart(&bp[k], props); + /* And the tracers */ + tracers_first_init_bpart(&bp[k], us, phys_const, cosmo); + /* And the splitting data */ particle_splitting_mark_part_as_not_split(&bp[k].split_data, bp[k].id); diff --git a/src/stars/Basic/stars_part.h b/src/stars/Basic/stars_part.h index 4fa337e2020e25f63c9db8e8dc0eb866e422b343..713fe0e281c65ad4e484e42682f83c16b88b1cda 100644 --- a/src/stars/Basic/stars_part.h +++ b/src/stars/Basic/stars_part.h @@ -88,7 +88,7 @@ struct spart { struct feedback_spart_data feedback_data; /*! Tracer structure */ - struct tracers_xpart_data tracers_data; + struct tracers_spart_data tracers_data; /*! Chemistry structure */ struct chemistry_spart_data chemistry_data; diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h index b00abae2fe0ac69fdfc8b6a153a842b9b7bb3227..22091ecc4f6c43928993b6eb8ff43ed9d9e7cd1e 100644 --- a/src/stars/EAGLE/stars_part.h +++ b/src/stars/EAGLE/stars_part.h @@ -110,7 +110,7 @@ struct spart { struct feedback_spart_data feedback_data; /*! Tracer structure */ - struct tracers_xpart_data tracers_data; + struct tracers_spart_data tracers_data; /*! Chemistry structure */ struct chemistry_spart_data chemistry_data; diff --git a/src/stars/GEAR/stars_part.h b/src/stars/GEAR/stars_part.h index 0e4406110bf09dfaee15f0fd815746901169a06c..d9d73c1f3e2da40aeef359e576a3d9b298b10a5b 100644 --- a/src/stars/GEAR/stars_part.h +++ b/src/stars/GEAR/stars_part.h @@ -88,7 +88,7 @@ struct spart { struct feedback_spart_data feedback_data; /*! Tracer structure */ - struct tracers_xpart_data tracers_data; + struct tracers_spart_data tracers_data; /*! Chemistry structure */ struct chemistry_spart_data chemistry_data; diff --git a/src/stars/None/stars_part.h b/src/stars/None/stars_part.h index 94fcf5597de6a221bf955faa1d200b31687ae116..2400407ce25a60b504c31bb782730a7dfe6521fc 100644 --- a/src/stars/None/stars_part.h +++ b/src/stars/None/stars_part.h @@ -58,6 +58,9 @@ struct spart { /*! Particle time bin */ timebin_t time_bin; + /*! Tracer structure */ + struct tracers_spart_data tracers_data; + /*! Splitting structure */ struct particle_splitting_data split_data; diff --git a/src/tracers.h b/src/tracers.h index 9e939207c9ee98818fc9e788f05f5ae7662a961f..98ee03111561aeb8f559108c1d22c142425b9366 100644 --- a/src/tracers.h +++ b/src/tracers.h @@ -27,7 +27,7 @@ /* Config parameters. */ #include <config.h> -/* Import the right cooling definition */ +/* Import the right tracers definition */ #if defined(TRACERS_NONE) #include "./tracers/none/tracers.h" #elif defined(TRACERS_EAGLE) diff --git a/src/tracers/EAGLE/tracers.h b/src/tracers/EAGLE/tracers.h index 6ca4f1adfb62c2bf0d530e16ef22f83e18dfb711..6f135882acdcc1077c645228d154e9085ffcec74 100644 --- a/src/tracers/EAGLE/tracers.h +++ b/src/tracers/EAGLE/tracers.h @@ -24,8 +24,11 @@ #include <config.h> /* Local includes */ +#include "black_holes.h" #include "cooling.h" +#include "engine.h" #include "part.h" +#include "star_formation.h" #include "tracers_struct.h" /** @@ -77,7 +80,7 @@ static INLINE void tracers_after_drift( * @brief Update the particle tracers just after its time-step has been * computed. * - * In EAGLE we record the highest temperature reached. + * In EAGLE we record the highest temperature reached and the average SFR. * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data (containing the tracers @@ -89,12 +92,16 @@ static INLINE void tracers_after_drift( * @param hydro_props the hydro_props struct * @param cooling The #cooling_function_data used in the run. * @param time The current time. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_part) */ -static INLINE void tracers_after_timestep( +static INLINE void tracers_after_timestep_part( const struct part *p, struct xpart *xp, const struct unit_system *us, const struct phys_const *phys_const, const int with_cosmology, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct cooling_function_data *cooling, const double time) { + const struct cooling_function_data *cooling, const double time, + const double time_step_length, const int *const tracers_triggers_started) { /* Current temperature */ const float temperature = cooling_get_temperature(phys_const, hydro_props, us, @@ -111,6 +118,69 @@ static INLINE void tracers_after_timestep( xp->tracers_data.maximum_temperature_time = time; } } + + /* Accumulate average SFR */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + if (tracers_triggers_started[i]) + xp->tracers_data.averaged_SFR[i] += + star_formation_get_SFR(p, xp) * time_step_length; + } +} + +/** + * @brief Update the star particle tracers just after its time-step has been + * computed. + * + * In EAGLE, nothing to do. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param with_cosmology Are we running a cosmological simulation? + * @param cosmo The current cosmological model. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_spart) + */ +static INLINE void tracers_after_timestep_spart( + struct spart *sp, const struct unit_system *us, + const struct phys_const *phys_const, const int with_cosmology, + const struct cosmology *cosmo, const double time_step_length, + const int *const tracers_triggers_started) {} + +/** + * @brief Update the black hole particle tracers just after its time-step has + * been computed. + * + * In EAGLE, we record the average accr. rate. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param with_cosmology Are we running a cosmological simulation? + * @param cosmo The current cosmological model. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_bpart) + */ +static INLINE void tracers_after_timestep_bpart( + struct bpart *bp, const struct unit_system *us, + const struct phys_const *phys_const, const int with_cosmology, + const struct cosmology *cosmo, const double time_step_length, + const int *const tracers_triggers_started) { + + const float accr_rate = black_holes_get_accretion_rate(bp); + + /* Accumulate average accretion rate */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + if (tracers_triggers_started[i]) + bp->tracers_data.averaged_accretion_rate[i] += + accr_rate * time_step_length; + } } /** @@ -146,6 +216,41 @@ static INLINE void tracers_first_init_xpart( xp->tracers_data.density_at_last_AGN_feedback_event = -1.f; } +/** + * @brief Initialise the star tracer data at the start of a calculation. + * + * In EAGLE, nothing to do. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + */ +static INLINE void tracers_first_init_spart(struct spart *sp, + const struct unit_system *us, + const struct phys_const *phys_const, + const struct cosmology *cosmo) {} + +/** + * @brief Initialise the black hole tracer data at the start of a calculation. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + */ +static INLINE void tracers_first_init_bpart(struct bpart *bp, + const struct unit_system *us, + const struct phys_const *phys_const, + const struct cosmology *cosmo) { + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) + bp->tracers_data.averaged_accretion_rate[i] = 0.f; +} + /** * @brief Update the particles' tracer data after a stellar feedback * event. @@ -222,6 +327,41 @@ static INLINE void tracers_after_jet_feedback( xp->tracers_data.jet_feedback_energy += delta_energy; } +/** + * @brief Tracer event called after a snapshot was written. + * + * @param p the #part. + * @param xp the #xpart. + */ +static INLINE void tracers_after_snapshot_part(const struct part *p, + struct xpart *xp) { + + for (int i = 0; i < num_snapshot_triggers_part; ++i) + xp->tracers_data.averaged_SFR[i] = 0.f; +} + +/** + * @brief Tracer event called after a snapshot was written. + * + * @param sp the #spart. + */ +static INLINE void tracers_after_snapshot_spart(struct spart *sp) { + + for (int i = 0; i < num_snapshot_triggers_part; ++i) + sp->tracers_data.averaged_SFR[i] = 0.f; +} + +/** + * @brief Tracer event called after a snapshot was written. + * + * @param bp the #bpart. + */ +static INLINE void tracers_after_snapshot_bpart(struct bpart *bp) { + + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) + bp->tracers_data.averaged_accretion_rate[i] = 0.f; +} + /** * @brief Split the tracer content of a particle into n pieces * diff --git a/src/tracers/EAGLE/tracers_io.h b/src/tracers/EAGLE/tracers_io.h index a8b35458409e53bd9a041e003049bc86388a1f1a..9997449151cba56ffd542c2c2d9f28a5e1586e60 100644 --- a/src/tracers/EAGLE/tracers_io.h +++ b/src/tracers/EAGLE/tracers_io.h @@ -41,6 +41,48 @@ __attribute__((always_inline)) INLINE static void tracers_write_flavour( } #endif +INLINE static void convert_part_averaged_SFR(const struct engine* e, + const struct part* p, + const struct xpart* xp, + float* ret) { + + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + if (e->snapshot_recording_triggers_started_part[i]) + ret[i] = xp->tracers_data.averaged_SFR[i] / + e->snapshot_recording_triggers_part[i]; + else + ret[i] = 0.f; + } +} + +INLINE static void convert_spart_averaged_SFR(const struct engine* e, + const struct spart* sp, + float* ret) { + + /* Note: We use the 'part' trigger here as the SF would have started when the + * particle was still in gas form */ + for (int i = 0; i < num_snapshot_triggers_part; ++i) { + if (e->snapshot_recording_triggers_started_part[i]) + ret[i] = sp->tracers_data.averaged_SFR[i] / + e->snapshot_recording_triggers_part[i]; + else + ret[i] = 0.f; + } +} + +INLINE static void convert_bpart_averaged_accretion_rate(const struct engine* e, + const struct bpart* bp, + float* ret) { + + for (int i = 0; i < num_snapshot_triggers_bpart; ++i) { + if (e->snapshot_recording_triggers_started_bpart[i]) + ret[i] = bp->tracers_data.averaged_accretion_rate[i] / + e->snapshot_recording_triggers_bpart[i]; + else + ret[i] = 0.f; + } +} + /** * @brief Specifies which particle fields to write to a dataset * @@ -139,7 +181,13 @@ __attribute__((always_inline)) INLINE static int tracers_write_particles( "particle has never been hit by feedback"); } - return 10; + list[10] = io_make_output_field_convert_part( + "AveragedStarFormationRates", FLOAT, 2, UNIT_CONV_SFR, 0.f, parts, xparts, + convert_part_averaged_SFR, + "Star formation rates of the particles averaged over the period set by " + "the first two snapshot triggers"); + + return 11; } __attribute__((always_inline)) INLINE static int tracers_write_sparticles( @@ -237,7 +285,28 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles( "-1 if a particle has never been hit by feedback"); } - return 10; + list[10] = io_make_output_field_convert_spart( + "AveragedStarFormationRates", FLOAT, num_snapshot_triggers_part, + UNIT_CONV_SFR, 0.f, sparts, convert_spart_averaged_SFR, + "Star formation rates of the particles averaged over the period set by " + "the first two snapshot triggers when the particle was still a gas " + "particle."); + + return 11; +} + +__attribute__((always_inline)) INLINE static int tracers_write_bparticles( + const struct bpart* bparts, struct io_props* list, + const int with_cosmology) { + + list[0] = io_make_output_field_convert_bpart( + "AveragedAccretionRates", FLOAT, num_snapshot_triggers_bpart, + UNIT_CONV_MASS_PER_UNIT_TIME, 0.f, bparts, + convert_bpart_averaged_accretion_rate, + "Accretion rates of the black holes averaged over the period set by the " + "first two snapshot triggers"); + + return 1; } #endif /* SWIFT_TRACERS_EAGLE_IO_H */ diff --git a/src/tracers/EAGLE/tracers_struct.h b/src/tracers/EAGLE/tracers_struct.h index 71b1966a73c8515c725853a55dc2d87f12750a48..d1203ddd9fc54dd8d050c966da119873b64ab649 100644 --- a/src/tracers/EAGLE/tracers_struct.h +++ b/src/tracers/EAGLE/tracers_struct.h @@ -20,6 +20,9 @@ #ifndef SWIFT_TRACERS_STRUCT_EAGLE_H #define SWIFT_TRACERS_STRUCT_EAGLE_H +/* Local includes */ +#include "tracers_triggers.h" + /** * @brief Properties of the tracers stored in the extended particle data. */ @@ -56,6 +59,9 @@ struct tracers_xpart_data { float last_AGN_jet_feedback_time; }; + /*! Averaged SFR over two different time slices */ + float averaged_SFR[num_snapshot_triggers_part]; + /*! Density of the gas before the last AGN feedback event * (physical internal units) */ float density_before_last_AGN_feedback_event; @@ -90,4 +96,20 @@ struct tracers_xpart_data { char hit_by_AGN_feedback; }; +/** + * @brief Properties of the tracers stored in the star particle data. + * + * Note: In this model, they are identical to the xpart data. + */ +#define tracers_spart_data tracers_xpart_data + +/** + * @brief Properties of the tracers stored in the black hole particle data. + */ +struct tracers_bpart_data { + + /*! Averaged accretion rate over two different time slices */ + float averaged_accretion_rate[num_snapshot_triggers_bpart]; +}; + #endif /* SWIFT_TRACERS_STRUCT_EAGLE_H */ diff --git a/src/tracers/none/tracers.h b/src/tracers/none/tracers.h index 11c1028a15b39eb68006ced9fe1ec007d9fdf8a7..29c0017def96d2565c3e47c2394fde5c771ccc19 100644 --- a/src/tracers/none/tracers.h +++ b/src/tracers/none/tracers.h @@ -72,22 +72,74 @@ static INLINE void tracers_after_drift( * @brief Update the particle tracers just after its time-step has been * computed. * - * Nothing to do here in the EAGLE model. + * Nothing to do here. * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). * @param us The internal system of units. * @param phys_const The physical constants in internal units. + * @param with_cosmology Are we running a cosmological simulation? * @param cosmo The current cosmological model. * @param hydro_props the hydro_props struct * @param cooling The #cooling_function_data used in the run. + * @param time The current time. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_part) + */ +static INLINE void tracers_after_timestep_part( + const struct part *p, struct xpart *xp, const struct unit_system *us, + const struct phys_const *phys_const, const int with_cosmology, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct cooling_function_data *cooling, const double time, + const double time_step_length, const int *const tracers_triggers_started) {} + +/** + * @brief Update the star particle tracers just after its time-step has been + * computed. + * + * Nothing to do here. + * * @param p Pointer to the particle data. * @param xp Pointer to the extended particle data (containing the tracers * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param with_cosmology Are we running a cosmological simulation? + * @param cosmo The current cosmological model. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_spart) */ -static INLINE void tracers_after_timestep( - const struct part *p, struct xpart *xp, const struct unit_system *us, +static INLINE void tracers_after_timestep_spart( + struct spart *sp, const struct unit_system *us, const struct phys_const *phys_const, const int with_cosmology, - const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct cooling_function_data *cooling, const double time) {} + const struct cosmology *cosmo, const double time_step_length, + const int *const tracers_triggers_started) {} + +/** + * @brief Update the black hole particle tracers just after its time-step has + * been computed. + * + * Nothing to do here. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param with_cosmology Are we running a cosmological simulation? + * @param cosmo The current cosmological model. + * @param time_step_length The length of the step that just finished + * @param tracers_triggers_started Which triggers have started? (array of size + * num_snapshot_triggers_bpart) + */ +static INLINE void tracers_after_timestep_bpart( + struct bpart *bp, const struct unit_system *us, + const struct phys_const *phys_const, const int with_cosmology, + const struct cosmology *cosmo, const double time_step_length, + const int *const tracers_triggers_started) {} /** * @brief Initialise the tracer data at the start of a calculation. @@ -107,6 +159,40 @@ static INLINE void tracers_first_init_xpart( const struct hydro_props *hydro_props, const struct cooling_function_data *cooling) {} +/** + * @brief Initialise the star tracer data at the start of a calculation. + * + * Nothing to do here. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + */ +static INLINE void tracers_first_init_spart(struct spart *sp, + const struct unit_system *us, + const struct phys_const *phys_const, + const struct cosmology *cosmo) {} + +/** + * @brief Initialise the black hole tracer data at the start of a calculation. + * + * Nothing to do here. + * + * @param p Pointer to the particle data. + * @param xp Pointer to the extended particle data (containing the tracers + * struct). + * @param us The internal system of units. + * @param phys_const The physical constants in internal units. + * @param cosmo The current cosmological model. + */ +static INLINE void tracers_first_init_bpart(struct bpart *bp, + const struct unit_system *us, + const struct phys_const *phys_const, + const struct cosmology *cosmo) {} + /** * @brief Update the particles' tracer data after a stellar feedback * event. @@ -149,6 +235,35 @@ static INLINE void tracers_after_black_holes_feedback( const struct part *p, struct xpart *xp, const int with_cosmology, const float scale_factor, const double time, const double delta_energy) {} +/** + * @brief Tracer event called after a snapshot was written. + * + * Nothing to do here. + * + * @param p the #part. + * @param xp the #xpart. + */ +static INLINE void tracers_after_snapshot_part(const struct part *p, + struct xpart *xp) {} + +/** + * @brief Tracer event called after a snapshot was written. + * + * Nothing to do here. + * + * @param sp the #spart. + */ +static INLINE void tracers_after_snapshot_spart(struct spart *sp) {} + +/** + * @brief Tracer event called after a snapshot was written. + * + * Nothing to do here. + * + * @param sp the #spart. + */ +static INLINE void tracers_after_snapshot_bpart(struct bpart *bp) {} + /** * @brief Split the tracer content of a particle into n pieces * diff --git a/src/tracers/none/tracers_io.h b/src/tracers/none/tracers_io.h index 2421ff725b48c914a34cd79adf9f27fcf96b4301..a5fa95eaa7807b71af8f0333951dca8d562380d4 100644 --- a/src/tracers/none/tracers_io.h +++ b/src/tracers/none/tracers_io.h @@ -63,4 +63,12 @@ __attribute__((always_inline)) INLINE static int tracers_write_sparticles( return 0; } + +__attribute__((always_inline)) INLINE static int tracers_write_bparticles( + const struct bpart* bparts, struct io_props* list, + const int with_cosmology) { + + return 0; +} + #endif /* SWIFT_TRACERS_NONE_IO_H */ diff --git a/src/tracers/none/tracers_struct.h b/src/tracers/none/tracers_struct.h index 5ee0027aa13c76cf5ee73a128d36c34641c0030e..16c5d273a8062c3f31f6c35d1dba8f69fa3d1bcd 100644 --- a/src/tracers/none/tracers_struct.h +++ b/src/tracers/none/tracers_struct.h @@ -24,4 +24,14 @@ */ struct tracers_xpart_data {}; +/** + * @brief Properties of the tracers stored in the star particle data. + */ +struct tracers_spart_data {}; + +/** + * @brief Properties of the tracers stored in the black hole particle data. + */ +struct tracers_bpart_data {}; + #endif /* SWIFT_TRACERS_STRUCT_NONE_H */ diff --git a/src/tracers_triggers.h b/src/tracers_triggers.h new file mode 100644 index 0000000000000000000000000000000000000000..6d171208f720f3f30a6f4546f83d1a82bf43479c --- /dev/null +++ b/src/tracers_triggers.h @@ -0,0 +1,43 @@ +/******************************************************************************* + * This file is part of SWIFT. + * Copyright (c) 2022 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/>. + * + ******************************************************************************/ +#ifndef SWIFT_TRACERS_TRIGGERS_H +#define SWIFT_TRACERS_TRIGGERS_H + +/** + * @file src/tracers_triggers.h + * @brief Branches between the different particle data tracers + */ + +/* Config parameters. */ +#include <config.h> + +/* Import the right tracers definition */ +#if defined(TRACERS_NONE) +#define num_snapshot_triggers_part 0 +#define num_snapshot_triggers_spart 0 +#define num_snapshot_triggers_bpart 0 +#elif defined(TRACERS_EAGLE) +#define num_snapshot_triggers_part 2 +#define num_snapshot_triggers_spart 0 +#define num_snapshot_triggers_bpart 2 +#else +#error "Invalid choice of tracers." +#endif + +#endif /* SWIFT_TRACERS_TRIGGERS_H */