diff --git a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
index 4b807846de2871e1a6c9e7baf4793c9fd6475a9c..3270af86afffba13cfac6be7c7c21dd84980adc5 100644
--- a/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_ICs/EAGLE_12/eagle_12.yml
@@ -164,7 +164,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_ICs/EAGLE_12/vrconfig_6dfofbound_subhalos_SO_hydro.cfg b/examples/EAGLE_ICs/EAGLE_12/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..45d6a8e7b7ecfeba7832e2dc19b5c11411374c81
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_12/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
@@ -0,0 +1,186 @@
+#Configuration file for analysing DM particles (in either DM only simulation or in full Hydro)
+#runs 6DFOF + substructure algorithm, demands subhalos and FOF halos be self-bound, calculates many properties
+#Units currently set to take in as input, Mpc, 1e10 solar masses, km/s, output in same units
+#To set temporally unique halo ids, alter Snapshot_value=SNAP to appropriate value. Ie: for snapshot 12, change SNAP to 12
+
+################################
+#unit options, should always be provided
+################################
+#EDIT THIS SECTION!!!!
+
+#naming convention is EAGLE Hydro. This config will be generalised to allow specific naming conventions to be implemented besides predefined ones
+HDF_name_convention=6
+Input_includes_star_particle=1 #include star particles in hydro input
+Input_includes_bh_particle=1 #include bh particles in hydro input
+Input_includes_wind_particle=0 #include wind particles in hydro input (used by Illustris and moves particle type 0 to particle type 3 when decoupled from hydro forces). Here shown as example
+Input_includes_tracer_particle=0 #include tracer particles in hydro input (used by Illustris). Here shown as example
+Input_includes_extradm_particle=0 #include extra dm particles stored in particle type 2 and type 3, useful for zooms
+
+
+#units conversion from input input to desired internal unit
+Length_unit=1.0 #default code unit,
+Velocity_unit=1.0 #default velocity unit,
+Mass_unit=1.0 #default mass unit,
+#assumes input is in 1e10 msun, Mpc and km/s and output units are the same
+Gravity=43.0211349 #for 1e10 Msun, km/s and Mpc
+Hubble_unit=100.0 # assuming units are km/s and Mpc, then value of Hubble in km/s/Mpc
+
+#conversion of output length units to kpc
+Length_unit_to_kpc=1000.0
+#conversion of output velocity units to km/s
+Velocity_to_kms=1.0
+#conversion of output mass units to solar masses
+Mass_to_solarmass=1.0e10
+
+################################
+#input related
+################################
+#input is from a cosmological so can use parameters like box size, h, Omega_m to calculate length and density scales
+Cosmological_input=1
+#sets the total buffer size in bytes used to store temporary particle information
+#of mpi read threads before they are broadcast to the appropriate waiting non-read threads
+#if not set, default value is equivalent to 1e6 particles per mpi process, quite large
+#but significantly minimises the number of send/receives
+#in this example the buffer size is roughly that for a send/receive of 10000 particles
+#for 100 mpi processes
+MPI_particle_total_buf_size=100000000
+
+################################
+#search related options
+################################
+
+#how to search a simulation
+Particle_search_type=1 #search dark matter particles only
+#for baryon search
+Baryon_searchflag=2 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
+#for search for substruture
+Search_for_substructure=1 #if 0, end search once field objects are found
+#also useful for zoom simulations or simulations of individual objects, setting this flag means no field structure search is run
+Singlehalo_search=0 #if file is single halo in which one wishes to search for substructure. Here disabled.
+#additional option for field haloes
+Keep_FOF=0 #if field 6DFOF search is done, allows to keep structures found in 3DFOF (can be interpreted as the inter halo stellar mass when only stellar search is used).\n
+
+#minimum size for structures
+Minimum_size=20 #min particles
+Minimum_halo_size=32 #if field halos have different minimum sizes, otherwise set to -1.
+
+#for field fof halo search
+FoF_Field_search_type=3 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each initial 3D FOF candidate
+Halo_3D_linking_length=0.20
+
+#for mean field estimates and local velocity density distribution funciton estimator related quantiites, rarely need to change this
+Local_velocity_density_approximate_calculation=1 #calculates velocity density using approximative (and quicker) near neighbour search
+Cell_fraction = 0.01 #fraction of field fof halo used to determine mean velocity distribution function. Typical values are ~0.005-0.02
+Grid_type=1 #normal entropy based grid, shouldn't have to change
+Nsearch_velocity=32 #number of velocity neighbours used to calculate local velocity distribution function. Typial values are ~32
+Nsearch_physical=256 #numerof physical neighbours from which the nearest velocity neighbour set is based. Typical values are 128-512
+
+#for substructure search, rarely ever need to change this
+FoF_search_type=1 #default phase-space FOF search. Don't really need to change
+Iterative_searchflag=1 #iterative substructure search, for substructure find initial candidate substructures with smaller linking lengths then expand search region
+Outlier_threshold=2.5 #outlier threshold for a particle to be considered residing in substructure, that is how dynamically distinct a particle is. Typical values are >2
+Substructure_physical_linking_length=0.10
+Velocity_ratio=2.0 #ratio of speeds used in phase-space FOF
+Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical values are ~10-30
+Velocity_linking_length=0.20 #where scaled by structure dispersion
+Significance_level=1.0 #how significant a substructure is relative to Poisson noise. Values >= 1 are fine.
+
+#for iterative substructure search, rarely ever need to change this
+Iterative_threshold_factor=1.0 #change in threshold value when using iterative search. Here no increase in threshold if iterative or not
+Iterative_linking_length_factor=2.0 #increase in final linking final iterative substructure search
+Iterative_Vratio_factor=1.0 #change in Vratio when using iterative search. no change in vratio
+Iterative_ThetaOp_factor=1.0 #change in velocity opening angle. no change in velocity opening angle
+
+#for checking for halo merger remnants, which are defined as large, well separated phase-space density maxima
+Halo_core_search=2 # searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
+#if searching for cores, linking lengths. likely does not need to change much
+Use_adaptive_core_search=0 #calculate dispersions in configuration & vel space to determine linking lengths
+Use_phase_tensor_core_growth=2 #use full stepped phase-space tensor assignment
+Halo_core_ellx_fac=0.7 #how linking lengths are changed when searching for local 6DFOF cores,
+Halo_core_ellv_fac=2.0 #how velocity lengths based on dispersions are changed when searching for local 6DFOF cores
+Halo_core_ncellfac=0.005 #fraction of total halo particle number setting min size of a local 6DFOF core
+Halo_core_num_loops=8 #number of loops to iteratively search for cores
+Halo_core_loop_ellx_fac=0.75 #how much to change the configuration space linking per iteration
+Halo_core_loop_ellv_fac=1.0 #how much to change the velocity space linking per iteration
+Halo_core_loop_elln_fac=1.2 #how much to change the min number of particles per iteration
+Halo_core_phase_significance=2.0 #how significant a core must be in terms of dispersions (sigma) significance
+
+################################
+#Unbinding options (VELOCIraptor is able to accurately identify tidal debris so particles need not be bound to a structure)
+################################
+#unbinding related items
+Unbind_flag=1 #run unbinding
+#objects must have particles that meet the allowed kinetic to potential ratio AND also have some total fraction that are completely bound.
+Unbinding_type=0
+#run unbinding of field structures, aka halos. This is useful for sams and 6DFOF halos but may not be useful if interested in 3DFOF mass functions.
+Bound_halos=0
+#alpha factor used to determine whether particle is "bound" alaph*T+W<0. For standard subhalo catalogues use >0.9 but if interested in tidal debris 0.2-0.5
+Allowed_kinetic_potential_ratio=0.95
+Min_bound_mass_frac=0.65 #minimum bound mass fraction
+#don't keep background potential when unbinding, faster than recalculating
+Keep_background_potential=1
+Softening_length=0.0
+Frac_pot_ref=1.0
+Min_npot_ref=20
+Kinetic_reference_frame_type=0
+#extra options in new unbinding optimisation
+Unbinding_max_unbound_removal_fraction_per_iteration=0.5
+Unbinding_max_unbound_fraction=0.95
+Unbinding_max_unbound_fraction_allowed=0.025
+
+
+################################
+#Calculation of properties related options
+################################
+Virial_density=500 #user defined virial overdensity. Note that 200 rho_c, 200 rho_m and BN98 are already calculated.
+#when calculating properties, for field objects calculate inclusive masses
+Inclusive_halo_masses=3 #calculate inclusive masses for halos using full Spherical overdensity apertures
+#ensures that output is physical and not comoving distances per little h
+Comoving_units=0
+#calculate more (sub)halo properties (like angular momentum in spherical overdensity apertures, both inclusive and exclusive)
+Extensive_halo_properties_output=1
+#calculate aperture masses
+Calculate_aperture_quantities=1
+Number_of_apertures=5
+Aperture_values_in_kpc=5,10,30,50,100,
+Number_of_projected_apertures=5
+Projected_aperture_values_in_kpc=5,10,30,50,100,
+#calculate radial profiles
+Calculate_radial_profiles=1
+Number_of_radial_profile_bin_edges=20
+#default radial normalisation log rad bins, normed by R200crit, Integer flag of 0 is log bins and R200crit norm.
+Radial_profile_norm=0
+Radial_profile_bin_edges=-2.,-1.87379263,-1.74758526,-1.62137789,-1.49517052,-1.36896316,-1.24275579,-1.11654842,-0.99034105,-0.86413368,-0.73792631,-0.61171894,-0.48551157,-0.3593042,-0.23309684,-0.10688947,0.0193179,0.14552527,0.27173264,0.39794001,
+Iterate_cm_flag=0 #do not interate to determine centre-of-mass
+Sort_by_binding_energy=1 #sort particles by binding energy
+Reference_frame_for_properties=2 #use the minimum potential as reference frame about which to calculate properties
+
+################################
+#output related
+################################
+
+Write_group_array_file=0 #do not write a group array file
+Separate_output_files=0 #do not separate output into field and substructure files similar to subfind
+Binary_output=2 #Use HDF5 output (binary output 1, ascii 0, and HDF 2)
+#output particles residing in the spherical overdensity apertures of halos, only the particles exclusively belonging to halos
+Spherical_overdensity_halo_particle_list_output=1
+
+#halo ids are adjusted by this value * 1000000000000 (or 1000000 if code compiled with the LONGINTS option turned off)
+#to ensure that halo ids are temporally unique. So if you had 100 snapshots, for snap 100 set this to 100 and 100*1000000000000 will
+#be added to the halo id as set for this snapshot, so halo 1 becomes halo 100*1000000000000+1 and halo 1 of snap 0 would just have ID=1
+
+#ALTER THIS as part of a script to get temporally unique ids
+Snapshot_value=SNAP
+
+################################
+#other options
+################################
+Verbose=1 #how talkative do you want the code to be, 0 not much, 1 a lot, 2 chatterbox
+
+Metallicity_to_solarmetallicity=1.0 #conversion of output to solarmetallicity
+Star_formation_rate_to_solarmassperyear=1.0 #similar but for star formation rates
+Stellar_age_to_yr=1.0 #similar but for stellar ages
+Stellar_age_input_is_cosmological_scalefactor=1 #indicates stars store formation scalefactor
+Metallicity_input_unit_conversion_to_output_unit=1.0
+Stellar_age_input_unit_conversion_to_output_unit=1.0
+Star_formation_rate_input_unit_conversion_to_output_unit=1.0
diff --git a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
index 222e7a81ac7a979b6e9be9b4ac29c32342e80dae..6461982f70c9635652b8253887a9998777f37c69 100644
--- a/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_ICs/EAGLE_25/eagle_25.yml
@@ -165,7 +165,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_ICs/EAGLE_25/vrconfig_6dfofbound_subhalos_SO_hydro.cfg b/examples/EAGLE_ICs/EAGLE_25/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..45d6a8e7b7ecfeba7832e2dc19b5c11411374c81
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_25/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
@@ -0,0 +1,186 @@
+#Configuration file for analysing DM particles (in either DM only simulation or in full Hydro)
+#runs 6DFOF + substructure algorithm, demands subhalos and FOF halos be self-bound, calculates many properties
+#Units currently set to take in as input, Mpc, 1e10 solar masses, km/s, output in same units
+#To set temporally unique halo ids, alter Snapshot_value=SNAP to appropriate value. Ie: for snapshot 12, change SNAP to 12
+
+################################
+#unit options, should always be provided
+################################
+#EDIT THIS SECTION!!!!
+
+#naming convention is EAGLE Hydro. This config will be generalised to allow specific naming conventions to be implemented besides predefined ones
+HDF_name_convention=6
+Input_includes_star_particle=1 #include star particles in hydro input
+Input_includes_bh_particle=1 #include bh particles in hydro input
+Input_includes_wind_particle=0 #include wind particles in hydro input (used by Illustris and moves particle type 0 to particle type 3 when decoupled from hydro forces). Here shown as example
+Input_includes_tracer_particle=0 #include tracer particles in hydro input (used by Illustris). Here shown as example
+Input_includes_extradm_particle=0 #include extra dm particles stored in particle type 2 and type 3, useful for zooms
+
+
+#units conversion from input input to desired internal unit
+Length_unit=1.0 #default code unit,
+Velocity_unit=1.0 #default velocity unit,
+Mass_unit=1.0 #default mass unit,
+#assumes input is in 1e10 msun, Mpc and km/s and output units are the same
+Gravity=43.0211349 #for 1e10 Msun, km/s and Mpc
+Hubble_unit=100.0 # assuming units are km/s and Mpc, then value of Hubble in km/s/Mpc
+
+#conversion of output length units to kpc
+Length_unit_to_kpc=1000.0
+#conversion of output velocity units to km/s
+Velocity_to_kms=1.0
+#conversion of output mass units to solar masses
+Mass_to_solarmass=1.0e10
+
+################################
+#input related
+################################
+#input is from a cosmological so can use parameters like box size, h, Omega_m to calculate length and density scales
+Cosmological_input=1
+#sets the total buffer size in bytes used to store temporary particle information
+#of mpi read threads before they are broadcast to the appropriate waiting non-read threads
+#if not set, default value is equivalent to 1e6 particles per mpi process, quite large
+#but significantly minimises the number of send/receives
+#in this example the buffer size is roughly that for a send/receive of 10000 particles
+#for 100 mpi processes
+MPI_particle_total_buf_size=100000000
+
+################################
+#search related options
+################################
+
+#how to search a simulation
+Particle_search_type=1 #search dark matter particles only
+#for baryon search
+Baryon_searchflag=2 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
+#for search for substruture
+Search_for_substructure=1 #if 0, end search once field objects are found
+#also useful for zoom simulations or simulations of individual objects, setting this flag means no field structure search is run
+Singlehalo_search=0 #if file is single halo in which one wishes to search for substructure. Here disabled.
+#additional option for field haloes
+Keep_FOF=0 #if field 6DFOF search is done, allows to keep structures found in 3DFOF (can be interpreted as the inter halo stellar mass when only stellar search is used).\n
+
+#minimum size for structures
+Minimum_size=20 #min particles
+Minimum_halo_size=32 #if field halos have different minimum sizes, otherwise set to -1.
+
+#for field fof halo search
+FoF_Field_search_type=3 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each initial 3D FOF candidate
+Halo_3D_linking_length=0.20
+
+#for mean field estimates and local velocity density distribution funciton estimator related quantiites, rarely need to change this
+Local_velocity_density_approximate_calculation=1 #calculates velocity density using approximative (and quicker) near neighbour search
+Cell_fraction = 0.01 #fraction of field fof halo used to determine mean velocity distribution function. Typical values are ~0.005-0.02
+Grid_type=1 #normal entropy based grid, shouldn't have to change
+Nsearch_velocity=32 #number of velocity neighbours used to calculate local velocity distribution function. Typial values are ~32
+Nsearch_physical=256 #numerof physical neighbours from which the nearest velocity neighbour set is based. Typical values are 128-512
+
+#for substructure search, rarely ever need to change this
+FoF_search_type=1 #default phase-space FOF search. Don't really need to change
+Iterative_searchflag=1 #iterative substructure search, for substructure find initial candidate substructures with smaller linking lengths then expand search region
+Outlier_threshold=2.5 #outlier threshold for a particle to be considered residing in substructure, that is how dynamically distinct a particle is. Typical values are >2
+Substructure_physical_linking_length=0.10
+Velocity_ratio=2.0 #ratio of speeds used in phase-space FOF
+Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical values are ~10-30
+Velocity_linking_length=0.20 #where scaled by structure dispersion
+Significance_level=1.0 #how significant a substructure is relative to Poisson noise. Values >= 1 are fine.
+
+#for iterative substructure search, rarely ever need to change this
+Iterative_threshold_factor=1.0 #change in threshold value when using iterative search. Here no increase in threshold if iterative or not
+Iterative_linking_length_factor=2.0 #increase in final linking final iterative substructure search
+Iterative_Vratio_factor=1.0 #change in Vratio when using iterative search. no change in vratio
+Iterative_ThetaOp_factor=1.0 #change in velocity opening angle. no change in velocity opening angle
+
+#for checking for halo merger remnants, which are defined as large, well separated phase-space density maxima
+Halo_core_search=2 # searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
+#if searching for cores, linking lengths. likely does not need to change much
+Use_adaptive_core_search=0 #calculate dispersions in configuration & vel space to determine linking lengths
+Use_phase_tensor_core_growth=2 #use full stepped phase-space tensor assignment
+Halo_core_ellx_fac=0.7 #how linking lengths are changed when searching for local 6DFOF cores,
+Halo_core_ellv_fac=2.0 #how velocity lengths based on dispersions are changed when searching for local 6DFOF cores
+Halo_core_ncellfac=0.005 #fraction of total halo particle number setting min size of a local 6DFOF core
+Halo_core_num_loops=8 #number of loops to iteratively search for cores
+Halo_core_loop_ellx_fac=0.75 #how much to change the configuration space linking per iteration
+Halo_core_loop_ellv_fac=1.0 #how much to change the velocity space linking per iteration
+Halo_core_loop_elln_fac=1.2 #how much to change the min number of particles per iteration
+Halo_core_phase_significance=2.0 #how significant a core must be in terms of dispersions (sigma) significance
+
+################################
+#Unbinding options (VELOCIraptor is able to accurately identify tidal debris so particles need not be bound to a structure)
+################################
+#unbinding related items
+Unbind_flag=1 #run unbinding
+#objects must have particles that meet the allowed kinetic to potential ratio AND also have some total fraction that are completely bound.
+Unbinding_type=0
+#run unbinding of field structures, aka halos. This is useful for sams and 6DFOF halos but may not be useful if interested in 3DFOF mass functions.
+Bound_halos=0
+#alpha factor used to determine whether particle is "bound" alaph*T+W<0. For standard subhalo catalogues use >0.9 but if interested in tidal debris 0.2-0.5
+Allowed_kinetic_potential_ratio=0.95
+Min_bound_mass_frac=0.65 #minimum bound mass fraction
+#don't keep background potential when unbinding, faster than recalculating
+Keep_background_potential=1
+Softening_length=0.0
+Frac_pot_ref=1.0
+Min_npot_ref=20
+Kinetic_reference_frame_type=0
+#extra options in new unbinding optimisation
+Unbinding_max_unbound_removal_fraction_per_iteration=0.5
+Unbinding_max_unbound_fraction=0.95
+Unbinding_max_unbound_fraction_allowed=0.025
+
+
+################################
+#Calculation of properties related options
+################################
+Virial_density=500 #user defined virial overdensity. Note that 200 rho_c, 200 rho_m and BN98 are already calculated.
+#when calculating properties, for field objects calculate inclusive masses
+Inclusive_halo_masses=3 #calculate inclusive masses for halos using full Spherical overdensity apertures
+#ensures that output is physical and not comoving distances per little h
+Comoving_units=0
+#calculate more (sub)halo properties (like angular momentum in spherical overdensity apertures, both inclusive and exclusive)
+Extensive_halo_properties_output=1
+#calculate aperture masses
+Calculate_aperture_quantities=1
+Number_of_apertures=5
+Aperture_values_in_kpc=5,10,30,50,100,
+Number_of_projected_apertures=5
+Projected_aperture_values_in_kpc=5,10,30,50,100,
+#calculate radial profiles
+Calculate_radial_profiles=1
+Number_of_radial_profile_bin_edges=20
+#default radial normalisation log rad bins, normed by R200crit, Integer flag of 0 is log bins and R200crit norm.
+Radial_profile_norm=0
+Radial_profile_bin_edges=-2.,-1.87379263,-1.74758526,-1.62137789,-1.49517052,-1.36896316,-1.24275579,-1.11654842,-0.99034105,-0.86413368,-0.73792631,-0.61171894,-0.48551157,-0.3593042,-0.23309684,-0.10688947,0.0193179,0.14552527,0.27173264,0.39794001,
+Iterate_cm_flag=0 #do not interate to determine centre-of-mass
+Sort_by_binding_energy=1 #sort particles by binding energy
+Reference_frame_for_properties=2 #use the minimum potential as reference frame about which to calculate properties
+
+################################
+#output related
+################################
+
+Write_group_array_file=0 #do not write a group array file
+Separate_output_files=0 #do not separate output into field and substructure files similar to subfind
+Binary_output=2 #Use HDF5 output (binary output 1, ascii 0, and HDF 2)
+#output particles residing in the spherical overdensity apertures of halos, only the particles exclusively belonging to halos
+Spherical_overdensity_halo_particle_list_output=1
+
+#halo ids are adjusted by this value * 1000000000000 (or 1000000 if code compiled with the LONGINTS option turned off)
+#to ensure that halo ids are temporally unique. So if you had 100 snapshots, for snap 100 set this to 100 and 100*1000000000000 will
+#be added to the halo id as set for this snapshot, so halo 1 becomes halo 100*1000000000000+1 and halo 1 of snap 0 would just have ID=1
+
+#ALTER THIS as part of a script to get temporally unique ids
+Snapshot_value=SNAP
+
+################################
+#other options
+################################
+Verbose=1 #how talkative do you want the code to be, 0 not much, 1 a lot, 2 chatterbox
+
+Metallicity_to_solarmetallicity=1.0 #conversion of output to solarmetallicity
+Star_formation_rate_to_solarmassperyear=1.0 #similar but for star formation rates
+Stellar_age_to_yr=1.0 #similar but for stellar ages
+Stellar_age_input_is_cosmological_scalefactor=1 #indicates stars store formation scalefactor
+Metallicity_input_unit_conversion_to_output_unit=1.0
+Stellar_age_input_unit_conversion_to_output_unit=1.0
+Star_formation_rate_input_unit_conversion_to_output_unit=1.0
diff --git a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
index 27320d1327c412cc803eb3f8859a7f569bc14175..c6f7f0e01bd292b03b0eb31c6c1b57d5163dafe7 100644
--- a/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_ICs/EAGLE_50/eagle_50.yml
@@ -165,7 +165,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_ICs/EAGLE_50/vrconfig_6dfofbound_subhalos_SO_hydro.cfg b/examples/EAGLE_ICs/EAGLE_50/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..45d6a8e7b7ecfeba7832e2dc19b5c11411374c81
--- /dev/null
+++ b/examples/EAGLE_ICs/EAGLE_50/vrconfig_6dfofbound_subhalos_SO_hydro.cfg
@@ -0,0 +1,186 @@
+#Configuration file for analysing DM particles (in either DM only simulation or in full Hydro)
+#runs 6DFOF + substructure algorithm, demands subhalos and FOF halos be self-bound, calculates many properties
+#Units currently set to take in as input, Mpc, 1e10 solar masses, km/s, output in same units
+#To set temporally unique halo ids, alter Snapshot_value=SNAP to appropriate value. Ie: for snapshot 12, change SNAP to 12
+
+################################
+#unit options, should always be provided
+################################
+#EDIT THIS SECTION!!!!
+
+#naming convention is EAGLE Hydro. This config will be generalised to allow specific naming conventions to be implemented besides predefined ones
+HDF_name_convention=6
+Input_includes_star_particle=1 #include star particles in hydro input
+Input_includes_bh_particle=1 #include bh particles in hydro input
+Input_includes_wind_particle=0 #include wind particles in hydro input (used by Illustris and moves particle type 0 to particle type 3 when decoupled from hydro forces). Here shown as example
+Input_includes_tracer_particle=0 #include tracer particles in hydro input (used by Illustris). Here shown as example
+Input_includes_extradm_particle=0 #include extra dm particles stored in particle type 2 and type 3, useful for zooms
+
+
+#units conversion from input input to desired internal unit
+Length_unit=1.0 #default code unit,
+Velocity_unit=1.0 #default velocity unit,
+Mass_unit=1.0 #default mass unit,
+#assumes input is in 1e10 msun, Mpc and km/s and output units are the same
+Gravity=43.0211349 #for 1e10 Msun, km/s and Mpc
+Hubble_unit=100.0 # assuming units are km/s and Mpc, then value of Hubble in km/s/Mpc
+
+#conversion of output length units to kpc
+Length_unit_to_kpc=1000.0
+#conversion of output velocity units to km/s
+Velocity_to_kms=1.0
+#conversion of output mass units to solar masses
+Mass_to_solarmass=1.0e10
+
+################################
+#input related
+################################
+#input is from a cosmological so can use parameters like box size, h, Omega_m to calculate length and density scales
+Cosmological_input=1
+#sets the total buffer size in bytes used to store temporary particle information
+#of mpi read threads before they are broadcast to the appropriate waiting non-read threads
+#if not set, default value is equivalent to 1e6 particles per mpi process, quite large
+#but significantly minimises the number of send/receives
+#in this example the buffer size is roughly that for a send/receive of 10000 particles
+#for 100 mpi processes
+MPI_particle_total_buf_size=100000000
+
+################################
+#search related options
+################################
+
+#how to search a simulation
+Particle_search_type=1 #search dark matter particles only
+#for baryon search
+Baryon_searchflag=2 #if 1 search for baryons separately using phase-space search when identifying substructures, 2 allows special treatment in field FOF linking and phase-space substructure search, 0 treat the same as dark matter particles
+#for search for substruture
+Search_for_substructure=1 #if 0, end search once field objects are found
+#also useful for zoom simulations or simulations of individual objects, setting this flag means no field structure search is run
+Singlehalo_search=0 #if file is single halo in which one wishes to search for substructure. Here disabled.
+#additional option for field haloes
+Keep_FOF=0 #if field 6DFOF search is done, allows to keep structures found in 3DFOF (can be interpreted as the inter halo stellar mass when only stellar search is used).\n
+
+#minimum size for structures
+Minimum_size=20 #min particles
+Minimum_halo_size=32 #if field halos have different minimum sizes, otherwise set to -1.
+
+#for field fof halo search
+FoF_Field_search_type=3 #5 3DFOF search for field halos, 4 for 6DFOF clean up of field halos, 3 for 6DFOF with velocity scale distinct for each initial 3D FOF candidate
+Halo_3D_linking_length=0.20
+
+#for mean field estimates and local velocity density distribution funciton estimator related quantiites, rarely need to change this
+Local_velocity_density_approximate_calculation=1 #calculates velocity density using approximative (and quicker) near neighbour search
+Cell_fraction = 0.01 #fraction of field fof halo used to determine mean velocity distribution function. Typical values are ~0.005-0.02
+Grid_type=1 #normal entropy based grid, shouldn't have to change
+Nsearch_velocity=32 #number of velocity neighbours used to calculate local velocity distribution function. Typial values are ~32
+Nsearch_physical=256 #numerof physical neighbours from which the nearest velocity neighbour set is based. Typical values are 128-512
+
+#for substructure search, rarely ever need to change this
+FoF_search_type=1 #default phase-space FOF search. Don't really need to change
+Iterative_searchflag=1 #iterative substructure search, for substructure find initial candidate substructures with smaller linking lengths then expand search region
+Outlier_threshold=2.5 #outlier threshold for a particle to be considered residing in substructure, that is how dynamically distinct a particle is. Typical values are >2
+Substructure_physical_linking_length=0.10
+Velocity_ratio=2.0 #ratio of speeds used in phase-space FOF
+Velocity_opening_angle=0.10 #angle between velocities. 18 degrees here, typical values are ~10-30
+Velocity_linking_length=0.20 #where scaled by structure dispersion
+Significance_level=1.0 #how significant a substructure is relative to Poisson noise. Values >= 1 are fine.
+
+#for iterative substructure search, rarely ever need to change this
+Iterative_threshold_factor=1.0 #change in threshold value when using iterative search. Here no increase in threshold if iterative or not
+Iterative_linking_length_factor=2.0 #increase in final linking final iterative substructure search
+Iterative_Vratio_factor=1.0 #change in Vratio when using iterative search. no change in vratio
+Iterative_ThetaOp_factor=1.0 #change in velocity opening angle. no change in velocity opening angle
+
+#for checking for halo merger remnants, which are defined as large, well separated phase-space density maxima
+Halo_core_search=2 # searches for separate 6dfof cores in field haloes, and then more than just flags halo as merging, assigns particles to each merging "halo". 2 is full separation, 1 is flagging, 0 is off
+#if searching for cores, linking lengths. likely does not need to change much
+Use_adaptive_core_search=0 #calculate dispersions in configuration & vel space to determine linking lengths
+Use_phase_tensor_core_growth=2 #use full stepped phase-space tensor assignment
+Halo_core_ellx_fac=0.7 #how linking lengths are changed when searching for local 6DFOF cores,
+Halo_core_ellv_fac=2.0 #how velocity lengths based on dispersions are changed when searching for local 6DFOF cores
+Halo_core_ncellfac=0.005 #fraction of total halo particle number setting min size of a local 6DFOF core
+Halo_core_num_loops=8 #number of loops to iteratively search for cores
+Halo_core_loop_ellx_fac=0.75 #how much to change the configuration space linking per iteration
+Halo_core_loop_ellv_fac=1.0 #how much to change the velocity space linking per iteration
+Halo_core_loop_elln_fac=1.2 #how much to change the min number of particles per iteration
+Halo_core_phase_significance=2.0 #how significant a core must be in terms of dispersions (sigma) significance
+
+################################
+#Unbinding options (VELOCIraptor is able to accurately identify tidal debris so particles need not be bound to a structure)
+################################
+#unbinding related items
+Unbind_flag=1 #run unbinding
+#objects must have particles that meet the allowed kinetic to potential ratio AND also have some total fraction that are completely bound.
+Unbinding_type=0
+#run unbinding of field structures, aka halos. This is useful for sams and 6DFOF halos but may not be useful if interested in 3DFOF mass functions.
+Bound_halos=0
+#alpha factor used to determine whether particle is "bound" alaph*T+W<0. For standard subhalo catalogues use >0.9 but if interested in tidal debris 0.2-0.5
+Allowed_kinetic_potential_ratio=0.95
+Min_bound_mass_frac=0.65 #minimum bound mass fraction
+#don't keep background potential when unbinding, faster than recalculating
+Keep_background_potential=1
+Softening_length=0.0
+Frac_pot_ref=1.0
+Min_npot_ref=20
+Kinetic_reference_frame_type=0
+#extra options in new unbinding optimisation
+Unbinding_max_unbound_removal_fraction_per_iteration=0.5
+Unbinding_max_unbound_fraction=0.95
+Unbinding_max_unbound_fraction_allowed=0.025
+
+
+################################
+#Calculation of properties related options
+################################
+Virial_density=500 #user defined virial overdensity. Note that 200 rho_c, 200 rho_m and BN98 are already calculated.
+#when calculating properties, for field objects calculate inclusive masses
+Inclusive_halo_masses=3 #calculate inclusive masses for halos using full Spherical overdensity apertures
+#ensures that output is physical and not comoving distances per little h
+Comoving_units=0
+#calculate more (sub)halo properties (like angular momentum in spherical overdensity apertures, both inclusive and exclusive)
+Extensive_halo_properties_output=1
+#calculate aperture masses
+Calculate_aperture_quantities=1
+Number_of_apertures=5
+Aperture_values_in_kpc=5,10,30,50,100,
+Number_of_projected_apertures=5
+Projected_aperture_values_in_kpc=5,10,30,50,100,
+#calculate radial profiles
+Calculate_radial_profiles=1
+Number_of_radial_profile_bin_edges=20
+#default radial normalisation log rad bins, normed by R200crit, Integer flag of 0 is log bins and R200crit norm.
+Radial_profile_norm=0
+Radial_profile_bin_edges=-2.,-1.87379263,-1.74758526,-1.62137789,-1.49517052,-1.36896316,-1.24275579,-1.11654842,-0.99034105,-0.86413368,-0.73792631,-0.61171894,-0.48551157,-0.3593042,-0.23309684,-0.10688947,0.0193179,0.14552527,0.27173264,0.39794001,
+Iterate_cm_flag=0 #do not interate to determine centre-of-mass
+Sort_by_binding_energy=1 #sort particles by binding energy
+Reference_frame_for_properties=2 #use the minimum potential as reference frame about which to calculate properties
+
+################################
+#output related
+################################
+
+Write_group_array_file=0 #do not write a group array file
+Separate_output_files=0 #do not separate output into field and substructure files similar to subfind
+Binary_output=2 #Use HDF5 output (binary output 1, ascii 0, and HDF 2)
+#output particles residing in the spherical overdensity apertures of halos, only the particles exclusively belonging to halos
+Spherical_overdensity_halo_particle_list_output=1
+
+#halo ids are adjusted by this value * 1000000000000 (or 1000000 if code compiled with the LONGINTS option turned off)
+#to ensure that halo ids are temporally unique. So if you had 100 snapshots, for snap 100 set this to 100 and 100*1000000000000 will
+#be added to the halo id as set for this snapshot, so halo 1 becomes halo 100*1000000000000+1 and halo 1 of snap 0 would just have ID=1
+
+#ALTER THIS as part of a script to get temporally unique ids
+Snapshot_value=SNAP
+
+################################
+#other options
+################################
+Verbose=1 #how talkative do you want the code to be, 0 not much, 1 a lot, 2 chatterbox
+
+Metallicity_to_solarmetallicity=1.0 #conversion of output to solarmetallicity
+Star_formation_rate_to_solarmassperyear=1.0 #similar but for star formation rates
+Stellar_age_to_yr=1.0 #similar but for stellar ages
+Stellar_age_input_is_cosmological_scalefactor=1 #indicates stars store formation scalefactor
+Metallicity_input_unit_conversion_to_output_unit=1.0
+Stellar_age_input_unit_conversion_to_output_unit=1.0
+Star_formation_rate_input_unit_conversion_to_output_unit=1.0
diff --git a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
index faeda2a3acf4f04cd75512190666374351592276..9205b9f6f560398b046625f6c65fcc887492991d 100644
--- a/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_low_z/EAGLE_12/eagle_12.yml
@@ -156,7 +156,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
index 5736fb864bb950703c6861cefec063ceec7861ec..ccd89efbead0c511316bd664f6754ed0797c7c50 100644
--- a/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_low_z/EAGLE_25/eagle_25.yml
@@ -163,7 +163,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
index cb2fa3b4dc0b2e3f04f04f73b6157e1400665c67..e2a98725a94c9fd34b45a95530f480562c0aaec8 100644
--- a/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_low_z/EAGLE_50/eagle_50.yml
@@ -155,7 +155,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
index d1bbfa4ef6c6b785265bfda7a9d929374af5f962..9136f85424beaf70b614cde0e70b67c134056831 100644
--- a/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_low_z/EAGLE_6/eagle_6.yml
@@ -165,7 +165,7 @@ EAGLEFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
- max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ max_eddington_fraction: 1.5 # Maximal allowed accretion rate in units of the Eddington rate.
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
diff --git a/examples/SubgridTests/BlackHoleSwallowing/swallowing.yml b/examples/SubgridTests/BlackHoleSwallowing/swallowing.yml
index e0637308e9cc9631e7cd042920163a98bf7e0bf0..2390918132e73de031b3795d17d9a2cbe9863680 100644
--- a/examples/SubgridTests/BlackHoleSwallowing/swallowing.yml
+++ b/examples/SubgridTests/BlackHoleSwallowing/swallowing.yml
@@ -86,6 +86,7 @@ EAGLECooling:
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term.
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
diff --git a/examples/SubgridTests/BlackHoleSwallowing/swallowing_mpi.yml b/examples/SubgridTests/BlackHoleSwallowing/swallowing_mpi.yml
index 8db524c9caa7f0941e34d6a44976076abc7f2dba..c587a5b6ff9c125e51180a84f22844b0102bcb64 100644
--- a/examples/SubgridTests/BlackHoleSwallowing/swallowing_mpi.yml
+++ b/examples/SubgridTests/BlackHoleSwallowing/swallowing_mpi.yml
@@ -86,6 +86,7 @@ EAGLECooling:
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
+ viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
diff --git a/src/black_holes/Default/black_holes_io.h b/src/black_holes/Default/black_holes_io.h
index e193f550bfd100077c33abdf6bdfcedb74d829da..41ca7c1cd0b3542959d8a7372179c7e0cb285285 100644
--- a/src/black_holes/Default/black_holes_io.h
+++ b/src/black_holes/Default/black_holes_io.h
@@ -22,6 +22,20 @@
#include "black_holes_part.h"
#include "io_properties.h"
+INLINE static void convert_bpart_pos(const struct engine *e,
+ const struct bpart *bp, double *ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(bp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(bp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(bp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = bp->x[0];
+ ret[1] = bp->x[1];
+ ret[2] = bp->x[2];
+ }
+}
+
/**
* @brief Specifies which b-particle fields to read from a dataset
*
@@ -37,8 +51,9 @@ INLINE static void black_holes_read_particles(struct bpart *bparts,
*num_fields = 5;
/* List what we want to read */
- list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
- UNIT_CONV_LENGTH, bparts, x);
+ list[0] = io_make_output_field_convert_bpart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, bparts, convert_bpart_pos);
+
list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
UNIT_CONV_SPEED, bparts, v);
list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
diff --git a/src/black_holes/EAGLE/black_holes_io.h b/src/black_holes/EAGLE/black_holes_io.h
index cc6cb962db61f9e45ded990b9dca7f00db075733..f5a8da344ffe0d37564710b9ce9bc3e0d8f3d705 100644
--- a/src/black_holes/EAGLE/black_holes_io.h
+++ b/src/black_holes/EAGLE/black_holes_io.h
@@ -22,6 +22,20 @@
#include "black_holes_part.h"
#include "io_properties.h"
+INLINE static void convert_bpart_pos(const struct engine *e,
+ const struct bpart *bp, double *ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(bp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(bp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(bp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = bp->x[0];
+ ret[1] = bp->x[1];
+ ret[2] = bp->x[2];
+ }
+}
+
/**
* @brief Specifies which b-particle fields to read from a dataset
*
@@ -66,8 +80,8 @@ INLINE static void black_holes_write_particles(const struct bpart *bparts,
*num_fields = 12;
/* List what we want to write */
- list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
- bparts, x);
+ list[0] = io_make_output_field_convert_bpart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, bparts, convert_bpart_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, bparts, v);
list[2] =
diff --git a/src/cell.c b/src/cell.c
index 763f5117f17eec3927446e6a2dddab4b9b4ad4ec..3eb8b072de3fc26aefdfc6f4f8757ec10ee8280f 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -4324,23 +4324,34 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
/* In non-periodic BC runs, remove particles that crossed the border */
if (!periodic) {
+
/* Did the particle leave the box? */
if ((p->x[0] > dim[0]) || (p->x[0] < 0.) || // x
(p->x[1] > dim[1]) || (p->x[1] < 0.) || // y
(p->x[2] > dim[2]) || (p->x[2] < 0.)) { // z
- /* One last action before death? */
- hydro_remove_part(p, xp);
+ lock_lock(&e->s->lock);
- /* Remove the particle entirely */
- struct gpart *gp = p->gpart;
- cell_remove_part(e, c, p, xp);
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!part_is_inhibited(p, e)) {
- /* and it's gravity friend */
- if (gp != NULL) {
- cell_remove_gpart(e, c, gp);
+ /* One last action before death? */
+ hydro_remove_part(p, xp);
+
+ /* Remove the particle entirely */
+ struct gpart *gp = p->gpart;
+ cell_remove_part(e, c, p, xp);
+
+ /* and it's gravity friend */
+ if (gp != NULL) {
+ cell_remove_gpart(e, c, gp);
+ }
}
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
+
continue;
}
}
@@ -4484,13 +4495,24 @@ void cell_drift_gpart(struct cell *c, const struct engine *e, int force) {
/* In non-periodic BC runs, remove particles that crossed the border */
if (!periodic) {
+
/* Did the particle leave the box? */
if ((gp->x[0] > dim[0]) || (gp->x[0] < 0.) || // x
(gp->x[1] > dim[1]) || (gp->x[1] < 0.) || // y
(gp->x[2] > dim[2]) || (gp->x[2] < 0.)) { // z
- /* Remove the particle entirely */
- cell_remove_gpart(e, c, gp);
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!gpart_is_inhibited(gp, e)) {
+
+ /* Remove the particle entirely */
+ cell_remove_gpart(e, c, gp);
+ }
+
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
continue;
}
@@ -4614,17 +4636,28 @@ void cell_drift_spart(struct cell *c, const struct engine *e, int force) {
/* In non-periodic BC runs, remove particles that crossed the border */
if (!periodic) {
+
/* Did the particle leave the box? */
if ((sp->x[0] > dim[0]) || (sp->x[0] < 0.) || // x
(sp->x[1] > dim[1]) || (sp->x[1] < 0.) || // y
(sp->x[2] > dim[2]) || (sp->x[2] < 0.)) { // z
- /* Remove the particle entirely */
- struct gpart *gp = sp->gpart;
- cell_remove_spart(e, c, sp);
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!spart_is_inhibited(sp, e)) {
+
+ /* Remove the particle entirely */
+ struct gpart *gp = sp->gpart;
+ cell_remove_spart(e, c, sp);
+
+ /* and it's gravity friend */
+ cell_remove_gpart(e, c, gp);
+ }
- /* and it's gravity friend */
- cell_remove_gpart(e, c, gp);
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
continue;
}
@@ -4784,12 +4817,22 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
(bp->x[1] > dim[1]) || (bp->x[1] < 0.) || // y
(bp->x[2] > dim[2]) || (bp->x[2] < 0.)) { // z
- /* Remove the particle entirely */
- struct gpart *gp = bp->gpart;
- cell_remove_bpart(e, c, bp);
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!bpart_is_inhibited(bp, e)) {
+
+ /* Remove the particle entirely */
+ struct gpart *gp = bp->gpart;
+ cell_remove_bpart(e, c, bp);
+
+ /* and it's gravity friend */
+ cell_remove_gpart(e, c, gp);
+ }
- /* and it's gravity friend */
- cell_remove_gpart(e, c, gp);
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
continue;
}
diff --git a/src/part.c b/src/part.c
index 9394ab314b6c67cea832b2b5b96a8df76863ec8c..7d967caad18caf0c45d2a422fb4766d34f559d00 100644
--- a/src/part.c
+++ b/src/part.c
@@ -30,6 +30,7 @@
/* Local headers */
#include "error.h"
+#include "hydro.h"
#include "part.h"
/**
@@ -210,6 +211,13 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
part->x[1], part->x[2], gparts[k].x[0] - part->x[0],
gparts[k].x[1] - part->x[1], gparts[k].x[2] - part->x[2]);
+ /* Check that the particles have the same mass */
+ if (gparts[k].mass != hydro_get_mass(part))
+ error(
+ "Linked particles do not have the same mass!\n"
+ "gp->m=%e p->m=%e",
+ gparts[k].mass, hydro_get_mass(part));
+
/* Check that the particles are at the same time */
if (gparts[k].time_bin != part->time_bin)
error("Linked particles are not at the same time !");
@@ -237,6 +245,13 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
spart->x[1], spart->x[2], gparts[k].x[0] - spart->x[0],
gparts[k].x[1] - spart->x[1], gparts[k].x[2] - spart->x[2]);
+ /* Check that the particles have the same mass */
+ if (gparts[k].mass != spart->mass)
+ error(
+ "Linked particles do not have the same mass!\n"
+ "gp->m=%e sp->m=%e",
+ gparts[k].mass, spart->mass);
+
/* Check that the particles are at the same time */
if (gparts[k].time_bin != spart->time_bin)
error("Linked particles are not at the same time !");
@@ -264,6 +279,13 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
bpart->x[1], bpart->x[2], gparts[k].x[0] - bpart->x[0],
gparts[k].x[1] - bpart->x[1], gparts[k].x[2] - bpart->x[2]);
+ /* Check that the particles have the same mass */
+ if (gparts[k].mass != bpart->mass)
+ error(
+ "Linked particles do not have the same mass!\n"
+ "gp->m=%e sp->m=%e",
+ gparts[k].mass, bpart->mass);
+
/* Check that the particles are at the same time */
if (gparts[k].time_bin != bpart->time_bin)
error("Linked particles are not at the same time !");
@@ -287,6 +309,10 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
parts[k].x[2] != parts[k].gpart->x[2])
error("Linked particles are not at the same position !");
+ /* Check that the particles have the same mass */
+ if (hydro_get_mass(&parts[k]) != parts[k].gpart->mass)
+ error("Linked particles do not have the same mass!\n");
+
/* Check that the particles are at the same time */
if (parts[k].time_bin != parts[k].gpart->time_bin)
error("Linked particles are not at the same time !");
@@ -309,6 +335,10 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
sparts[k].x[2] != sparts[k].gpart->x[2])
error("Linked particles are not at the same position !");
+ /* Check that the particles have the same mass */
+ if (sparts[k].mass != sparts[k].gpart->mass)
+ error("Linked particles do not have the same mass!\n");
+
/* Check that the particles are at the same time */
if (sparts[k].time_bin != sparts[k].gpart->time_bin)
error("Linked particles are not at the same time !");
@@ -332,6 +362,10 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
bparts[k].x[2] != bparts[k].gpart->x[2])
error("Linked particles are not at the same position !");
+ /* Check that the particles have the same mass */
+ if (bparts[k].mass != bparts[k].gpart->mass)
+ error("Linked particles do not have the same mass!\n");
+
/* Check that the particles are at the same time */
if (bparts[k].time_bin != bparts[k].gpart->time_bin)
error("Linked particles are not at the same time !");
diff --git a/src/runner.c b/src/runner.c
index a196fa389a8d2b2edc290c8469dca503f8b2149d..bafb5fbf255a3c1c437eef04c090b8edd0522ddf 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -3819,10 +3819,20 @@ void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
message("BH %lld removing particle %lld", bp->id, p->id);
- /* Finally, remove the gas particle from the system */
- struct gpart *gp = p->gpart;
- cell_remove_part(e, c, p, xp);
- cell_remove_gpart(e, c, gp);
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!part_is_inhibited(p, e)) {
+
+ /* Finally, remove the gas particle from the system */
+ struct gpart *gp = p->gpart;
+ cell_remove_part(e, c, p, xp);
+ cell_remove_gpart(e, c, gp);
+ }
+
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
}
/* In any case, prevent the particle from being re-swallowed */
@@ -3852,10 +3862,20 @@ void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
message("BH %lld removing particle %lld (foreign BH case)",
bp->id, p->id);
- /* Finally, remove the gas particle from the system */
- struct gpart *gp = p->gpart;
- cell_remove_part(e, c, p, xp);
- cell_remove_gpart(e, c, gp);
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!part_is_inhibited(p, e)) {
+
+ /* Finally, remove the gas particle from the system */
+ struct gpart *gp = p->gpart;
+ cell_remove_part(e, c, p, xp);
+ cell_remove_gpart(e, c, gp);
+ }
+
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
found = 1;
break;
diff --git a/src/space.c b/src/space.c
index f622fda5f6e47205aeed2ce1f108848a81d2e00c..9b30d5cf176bc9f7ab7e64b331eb671a06e6b97e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -3936,7 +3936,7 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
for (int k = 0; k < nr_gparts; k++) {
/* Get the particle */
- const struct gpart *restrict gp = &gparts[k];
+ struct gpart *restrict gp = &gparts[k];
if (gp->type == swift_type_dark_matter)
continue;
@@ -3955,6 +3955,8 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
xp->v_full[0] = gp->v_full[0];
xp->v_full[1] = gp->v_full[1];
xp->v_full[2] = gp->v_full[2];
+
+ gp->mass = hydro_get_mass(p);
}
else if (gp->type == swift_type_stars) {
@@ -3966,6 +3968,8 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
sp->x[0] = gp->x[0];
sp->x[1] = gp->x[1];
sp->x[2] = gp->x[2];
+
+ gp->mass = sp->mass;
}
else if (gp->type == swift_type_black_hole) {
@@ -3977,6 +3981,8 @@ void space_synchronize_particle_positions_mapper(void *map_data, int nr_gparts,
bp->x[0] = gp->x[0];
bp->x[1] = gp->x[1];
bp->x[2] = gp->x[2];
+
+ gp->mass = bp->mass;
}
}
}
diff --git a/src/stars/Default/stars_io.h b/src/stars/Default/stars_io.h
index 5ff57549b9d0db37d6929c1e9fbb8da8372d7e6c..a8ec1cfa55728f9ca8a348d8fd6ec07d06b72185 100644
--- a/src/stars/Default/stars_io.h
+++ b/src/stars/Default/stars_io.h
@@ -22,6 +22,20 @@
#include "io_properties.h"
#include "stars_part.h"
+INLINE static void convert_spart_pos(const struct engine *e,
+ const struct spart *sp, double *ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(sp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(sp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(sp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = sp->x[0];
+ ret[1] = sp->x[1];
+ ret[2] = sp->x[2];
+ }
+}
+
/**
* @brief Specifies which s-particle fields to read from a dataset
*
@@ -64,8 +78,8 @@ INLINE static void stars_write_particles(const struct spart *sparts,
*num_fields = 5;
/* List what we want to write */
- list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
- sparts, x);
+ list[0] = io_make_output_field_convert_spart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, sparts, convert_spart_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
list[2] =
diff --git a/src/stars/EAGLE/stars_io.h b/src/stars/EAGLE/stars_io.h
index 0c03bee3007066c7c51c7ce0fb3d88d37a1b2ae3..cfacd52106398f435e56a9a2a67d1016726e2295 100644
--- a/src/stars/EAGLE/stars_io.h
+++ b/src/stars/EAGLE/stars_io.h
@@ -23,6 +23,20 @@
#include "io_properties.h"
#include "stars_part.h"
+INLINE static void convert_spart_pos(const struct engine *e,
+ const struct spart *sp, double *ret) {
+
+ if (e->s->periodic) {
+ ret[0] = box_wrap(sp->x[0], 0.0, e->s->dim[0]);
+ ret[1] = box_wrap(sp->x[1], 0.0, e->s->dim[1]);
+ ret[2] = box_wrap(sp->x[2], 0.0, e->s->dim[2]);
+ } else {
+ ret[0] = sp->x[0];
+ ret[1] = sp->x[1];
+ ret[2] = sp->x[2];
+ }
+}
+
/**
* @brief Specifies which s-particle fields to read from a dataset
*
@@ -67,8 +81,8 @@ INLINE static void stars_write_particles(const struct spart *sparts,
*num_fields = 10;
/* List what we want to write */
- list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
- sparts, x);
+ list[0] = io_make_output_field_convert_spart(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, sparts, convert_spart_pos);
list[1] =
io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, sparts, v);
list[2] =