...
 
Commits (78)
......@@ -185,7 +185,8 @@ EAGLEAGN:
max_eddington_fraction: 1.0 # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term
viscous_alpha: 1e6 # Normalisation constant of the Bondi viscuous time-scale accretion reduction term.
use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
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.
......
......@@ -209,6 +209,7 @@ EAGLEAGN:
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
......
......@@ -209,6 +209,7 @@ EAGLEAGN:
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
......
......@@ -207,6 +207,7 @@ EAGLEAGN:
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
......
......@@ -206,6 +206,7 @@ EAGLEAGN:
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
use_nibbling: 0 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
coupling_efficiency: 0.15 # Fraction of the radiated energy that couples to the gas in feedback events.
......
......@@ -18,6 +18,8 @@ Snapshots:
basename: pointMass # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
select_output_on: 1
select_output: select_output.yml
# Parameters governing the conserved quantities statistics
Statistics:
......
......@@ -532,28 +532,44 @@ GEARFeedback:
# EAGLE AGN model
EAGLEAGN:
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
use_subgrid_mass_from_ics: 1 # (Optional) Use subgrid masses specified in ICs [1, default], or initialise them to particle masses [0]?
with_subgrid_mass_check: 1 # (Optional) Verify that initial black hole subgrid masses are positive [1, default]. Only used if use_subgrid_mass_from_ics is 1.
multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH?
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
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.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
reposition_coefficient_upsilon: 0.0001 # Repositioning speed normalisation [km/s/M_sun]. Only meaningful if set_reposition_speed is 1.
reposition_exponent_xi: 1.0 # (Optional) Scaling of repositioning velocity with BH subgrid mass (default: 1.0, linear). Only meaningful if set_reposition_speed is 1.
threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
subgrid_seed_mass_Msun: 1.5e5 # Black hole subgrid mass at creation time in solar masses.
use_subgrid_mass_from_ics: 1 # (Optional) Use subgrid masses specified in ICs [1, default], or initialise them to particle masses [0]?
with_subgrid_mass_check: 1 # (Optional) Verify that initial black hole subgrid masses are positive [1, default]. Only used if use_subgrid_mass_from_ics is 1.
with_accretion_boost: 1 # Are we applying a Booth & Schaye (2009) like boost factor to black hole accretion rates?
accretion_boost_alpha: 1.0 # Constant pre-factor for BH accretion rates
accretion_boost_dens_norm: 0.1 # Normalisation density for boost factor [n_H / cm^3]; no boost is applied if the ambient density is lower. Only used if with_accretion_boost is 1.
accretion_boost_exponent: 2.0 # Exponent beta of boost factor model. Only used if with_accretion_boost is 1.
multi_phase_bondi: 0 # Compute Bondi rates per neighbour particle?
subgrid_bondi: 0 # Compute Bondi rates using the subgrid extrapolation of the gas properties around the BH?
with_angmom_limiter: 1 # Are we applying the Rosas-Guevara et al. (2015) viscous time-scale reduction term?
viscous_alpha: 1e6 # Normalisation constant of the viscous time-scale in the accretion reduction term
radiative_efficiency: 0.1 # Fraction of the accreted mass that gets radiated.
max_eddington_fraction: 1. # Maximal allowed accretion rate in units of the Eddington rate.
eddington_fraction_for_recording: 0.1 # Record the last time BHs reached an Eddington ratio above this threshold.
use_nibbling: 1 # Continuously transfer small amounts of mass from all gas neighbours to a black hole [1] or stochastically swallow whole gas particles [0]?
min_gas_mass_for_nibbling: 9e5 # Minimum mass for a gas particle to be nibbled from [M_Sun]. Only used if use_nibbling is 1.
use_scaled_coupling_efficiency: 1 # Switch to enable scaling of BH feedback efficiency with ambient gas density and/or metallicity.
coupling_efficiency: 0.15 # (Constant) fraction of the radiated energy that couples to the gas in feedback events. Only used if use_scaled_coupling_efficiency is 0.
epsilon_f_min: 0.01 # Minimal fraction of radiated black hole energy coupled to the gas. Only used if use_scaled_coupling_efficiency is 1.
epsilon_f_max: 0.2 # Maximal fraction of radiated black hole energy coupled to the gas. Only used if use_scaled_coupling_efficiency is 1.
epsilon_f_metallicity_norm: 0.0012663729 # Pivot point for the dependence of the black hole energy coupling efficiency on ambient gas metallicity (metal mass fraction). Only used if use_scaled_coupling_efficiency is 1.
epsilon_f_metallicity_exponent: 0.8686 # Power law index for the dependence of the black hole energy coupling efficiency on ambient gas metallicity. Only used if use_scaled_coupling_efficiency is 1.
epsilon_f_density_norm: 0.67 # Pivot point for the dependence of the black hole energy coupling efficiency on ambient gas density [n_H / cm^3]. Only used if use_scaled_coupling_efficiency is 1.
epsilon_f_density_exponent: 0.8686 # Power law index for the dependence of the black hole energy coupling efficiency on ambient gas density. Only used if use_scaled_coupling_efficiency is 1.
AGN_delta_T_K: 3.16228e8 # Change in temperature to apply to the gas particle in an AGN feedback event in Kelvin.
AGN_num_ngb_to_heat: 1. # Target number of gas neighbours to heat in an AGN feedback event.
max_reposition_mass: 2e8 # Maximal BH mass considered for BH repositioning in solar masses.
max_reposition_distance_ratio: 3.0 # Maximal distance a BH can be repositioned, in units of the softening length.
with_reposition_velocity_threshold: 1 # Should we only reposition to particles that move slowly w.r.t. the black hole?
max_reposition_velocity_ratio: 0.5 # Maximal velocity offset of a particle to reposition a BH to, in units of the ambient sound speed of the BH. Only meaningful if with_reposition_velocity_ratio is 1.
min_reposition_velocity_threshold: -1.0 # Minimal value of the velocity threshold for repositioning [km/s], set to < 0 for no effect. Only meaningful if with_reposition_velocity_ratio is 1.
set_reposition_speed: 0 # Should we reposition black holes with (at most) a prescribed speed towards the potential minimum?
reposition_coefficient_upsilon: 50.0 # Repositioning speed normalisation [km/s/M_sun]. Only meaningful if set_reposition_speed is 1.
reposition_reference_mass: 1e5 # Reference mass for scaling of reposition speed [M_sun]. Only meaningful if set_reposition_speed is 1.
reposition_exponent_mass: 2.0 # (Optional) Exponent for scaling of repositioning velocity with BH mass; default = 2.0. Only meaningful if set_reposition_speed is 1.
reposition_reference_n_H: 1.0 # Reference gas density around the black holes for scaling of repositioning speed [N_H cm^-3]. Only meaningful if set_reposition_speed is 1.
reposition_exponent_n_H: 1.0 # (Optional) Exponent for scaling of repositioning velocity with gas density; default = 1.0. Only meaningful if set_reposition_speed is 1.
threshold_major_merger: 0.333 # Mass ratio threshold to consider a BH merger as 'major'
threshold_minor_merger: 0.1 # Mass ratio threshold to consider a BH merger as 'minor'
merger_threshold_type: 2 # Type of velocity threshold for BH mergers (0: v_circ at kernel edge, 1: v_esc at actual distance, with softening, 2: v_esc at actual distance, no softening).
merger_max_distance_ratio: 3.0 # Maximal distance over which two BHs can merge, in units of the softening length.
......@@ -47,6 +47,10 @@ struct black_holes_props {
/*! Maximal change of h over one time-step */
float log_max_h_change;
/*! Use nibbling? */
int use_nibbling;
};
/**
......@@ -95,6 +99,8 @@ static INLINE void black_holes_props_init(struct black_holes_props *bp,
bp->log_max_h_change = hydro_props->log_max_h_change;
else
bp->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
bp->use_nibbling = 0;
}
/**
......
This diff is collapsed.
......@@ -81,11 +81,13 @@ runner_iact_nonsym_bh_gas_density(
/* Contribution to the total neighbour mass */
bi->ngb_mass += mj;
/* Neighbour sound speed */
const float cj = hydro_get_comoving_soundspeed(pj);
/* Contribution to the smoothed sound speed */
bi->sound_speed_gas += mj * cj * wi;
const double cj = hydro_get_comoving_soundspeed(pj);
bi->sound_speed_gas += mj * wi * cj;
/* Contribution to the smoothed gas metallicity */
bi->gas_metal_mass_fraction += mj * wi *
chemistry_get_metal_mass_fraction_for_black_holes(pj);
/* Neighbour internal energy */
const float uj = hydro_get_drifted_comoving_internal_energy(pj);
......@@ -106,9 +108,9 @@ runner_iact_nonsym_bh_gas_density(
/* Contribution to the specific angular momentum of gas, which is later
* converted to the circular velocity at the smoothing length */
bi->circular_velocity_gas[0] += mj * wi * (dx[1] * dv[2] - dx[2] * dv[1]);
bi->circular_velocity_gas[1] += mj * wi * (dx[2] * dv[0] - dx[0] * dv[2]);
bi->circular_velocity_gas[2] += mj * wi * (dx[0] * dv[1] - dx[1] * dv[0]);
bi->spec_angular_momentum_gas[0] += mj * wi * (dx[1] * dv[2] - dx[2] * dv[1]);
bi->spec_angular_momentum_gas[1] += mj * wi * (dx[2] * dv[0] - dx[0] * dv[2]);
bi->spec_angular_momentum_gas[2] += mj * wi * (dx[0] * dv[1] - dx[1] * dv[0]);
if (bh_props->multi_phase_bondi) {
/* Contribution to BH accretion rate
......@@ -256,8 +258,78 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
}
}
/* Is the BH hungry? */
if (bi->subgrid_mass > bi->mass) {
/* Check if the BH needs to be fed. If not, we're done here */
const float bh_mass_deficit = bi->subgrid_mass - bi->mass_at_start_of_step;
if (bh_mass_deficit <= 0) return;
if (bh_props->use_nibbling) {
/* If we do nibbling, things are quite straightforward. We transfer
* the mass and all associated quantities right here. */
if (bh_props->epsilon_r == 1)
return;
const float bi_mass_orig = bi->mass;
const float pj_mass_orig = hydro_get_mass(pj);
/* Don't nibble from particles that are too small already */
if (pj_mass_orig < bh_props->min_gas_mass_for_nibbling) return;
/* Next line is equivalent to w_ij * m_j / Sum_j (w_ij * m_j) */
const float particle_weight = hi_inv_dim * wi * pj_mass_orig / bi->rho_gas;
float nibble_mass = bh_mass_deficit * particle_weight;
/* We radiated away some of the accreted mass, so need to take slightly
* more from the gas than the BH gained */
const float excess_fraction = 1.0 / (1.0 - bh_props->epsilon_r);
/* Need to check whether nibbling would push gas mass below minimum
* allowed mass */
float new_gas_mass = pj_mass_orig - nibble_mass * excess_fraction;
if (new_gas_mass < bh_props->min_gas_mass_for_nibbling) {
new_gas_mass = bh_props->min_gas_mass_for_nibbling;
nibble_mass = (pj_mass_orig - bh_props->min_gas_mass_for_nibbling)
/ excess_fraction;
}
/* Transfer (dynamical) mass from the gas particle to the BH */
bi->mass += nibble_mass;
hydro_set_mass(pj, new_gas_mass);
/* Add the angular momentum of the accreted gas to the BH total.
* Note no change to gas here. The cosmological conversion factors for
* velocity (a^-1) and distance (a) cancel out, so the angular momentum
* is already in physical units. */
const float dv[3] = {bi->v[0] - pj->v[0], bi->v[1] - pj->v[1],
bi->v[2] - pj->v[2]};
bi->swallowed_angular_momentum[0] +=
nibble_mass * (dx[1] * dv[2] - dx[2] * dv[1]);
bi->swallowed_angular_momentum[1] +=
nibble_mass * (dx[2] * dv[0] - dx[0] * dv[2]);
bi->swallowed_angular_momentum[2] +=
nibble_mass * (dx[0] * dv[1] - dx[1] * dv[0]);
/* Update the BH momentum and velocity. Again, no change to gas here. */
const float bi_mom[3] = {bi_mass_orig * bi->v[0] + nibble_mass * pj->v[0],
bi_mass_orig * bi->v[1] + nibble_mass * pj->v[1],
bi_mass_orig * bi->v[2] + nibble_mass * pj->v[2]};
bi->v[0] = bi_mom[0] / bi->mass;
bi->v[1] = bi_mom[1] / bi->mass;
bi->v[2] = bi_mom[2] / bi->mass;
bi->gpart->v_full[0] = bi->v[0];
bi->gpart->v_full[1] = bi->v[1];
bi->gpart->v_full[2] = bi->v[2];
/* Update the BH and also gas metal masses */
struct chemistry_bpart_data *bi_chem = &bi->chemistry_data;
struct chemistry_part_data *pj_chem = &pj->chemistry_data;
chemistry_transfer_part_to_bpart(
bi_chem, pj_chem, nibble_mass * excess_fraction,
nibble_mass * excess_fraction / pj_mass_orig);
} else { /* ends nibbling section, below comes swallowing */
/* Probability to swallow this particle
* Recall that in SWIFT the SPH kernel is recovered by computing
......@@ -288,7 +360,7 @@ runner_iact_nonsym_bh_gas_swallow(const float r2, const float *dx,
bi->id, pj->id, pj->black_holes_data.swallow_id);
}
}
}
} /* ends section for swallowing */
}
/**
......
......@@ -21,6 +21,7 @@
#include "adiabatic_index.h"
#include "black_holes_part.h"
#include "black_holes_properties.h"
#include "io_properties.h"
/**
......@@ -125,6 +126,18 @@ INLINE static void convert_bpart_gas_circular_vel(const struct engine* e,
ret[2] = bp->circular_velocity_gas[2] * cosmo->a_inv;
}
INLINE static void convert_bpart_gas_temperatures(const struct engine* e,
const struct bpart* bp,
float* ret) {
const struct black_holes_props* props = e->black_holes_properties;
const struct cosmology* cosmo = e->cosmology;
/* Conversion from specific internal energy to temperature */
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
*
......@@ -139,7 +152,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
int with_cosmology) {
/* Say how much we want to write */
*num_fields = 30;
*num_fields = 44;
/* List what we want to write */
list[0] = io_make_output_field_convert_bpart(
......@@ -319,16 +332,110 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
"accretion disc model.");
list[27] = io_make_output_field(
"AccretedAngularMomenta", FLOAT, 3, UNIT_CONV_ANGULAR_MOMENTUM, 0.f,
bparts, accreted_angular_momentum,
"Physical angular momenta that the black holes have accumulated through "
"subgrid accretion.");
list[28] = io_make_output_field(
"CumulativeHeatingProbabilities", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
bparts, cumulative_target_prob,
"Cumulative (ideal) probability of the black holes heating any one of "
"their gas neighbour particles. This can be combined with "
"NumberOfTimeSteps to find the average heating probability between two "
"outputs.");
list[29] = io_make_output_field(
"CumulativeActualHeatingProbabilities", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
bparts, cumulative_actual_prob,
"Cumulative (actual) probability of the black holes heating any one of "
"their gas neighbour particles, accounting for the maximum allowed "
"probability. This can be combined with NumberOfTimeSteps to find the "
"average heating probability between two outputs.");
list[30] = io_make_output_field(
"HeatingProbabilities", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
bparts, target_heating_prob,
"Instantaneous (ideal) probability of the black holes heating any one of "
"their gas neighbour particles. The actual probability might be lower if "
"this exceeds the maximum allowed probability.");
list[31] = io_make_output_field(
"NumberOfGasNeighbours", INT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
num_ngbs,
"Integer number of gas neighbour particles within the black hole "
"kernels.");
list[32] = io_make_output_field(
"AccretionBoostFactors", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
accretion_boost_factor,
"Multiplicative factors by which the Bondi-Hoyle-Lyttleton accretion "
"rates have been increased by the density-dependent Booth & Schaye "
"(2009) accretion model.");
list[33] = io_make_output_field(
"GasMetallicities", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
gas_metal_mass_fraction,
"Metal mass fraction of the gas around the black holes.");
list[34] = io_make_output_field(
"FeedbackCouplingEfficiencies", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
epsilon_f,
"Fraction of available feedback energy that is coupled to the ambient "
"gas.");
list[35] = io_make_output_field(
"CumulativeFeedbackCouplingEfficiencies", FLOAT, 1, UNIT_CONV_NO_UNITS,
0.f, bparts, cumulative_epsilon_f,
"Cumulative fraction of available feedback energy that is coupled to "
"the ambient gas. This can be combined with NumberOfTimeSteps to find "
"the average coupling efficiency between two outputs.");
list[36] = io_make_output_field(
"FeedbackDeltaT", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, bparts,
AGN_delta_T,
"Temperature by which gas particles are heated by the black hole "
"particles in case of feedback.");
list[37] = io_make_output_field_convert_bpart(
"GasTemperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, bparts,
convert_bpart_gas_temperatures,
"Temperature of the gas surrounding the black holes.");
list[38] = io_make_output_field(
"SubgridDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts,
rho_subgrid_gas,
"Physical subgrid densities used in the subgrid-Bondi model.");
list[28] = io_make_output_field(
list[39] = io_make_output_field(
"SubgridSoundSpeeds", FLOAT, 1, UNIT_CONV_SPEED, 0.f, bparts,
sound_speed_subgrid_gas,
"Physical subgrid sound-speeds used in the subgrid-Bondi model.");
list[29] = io_make_output_field(
list[40] = io_make_output_field(
"LastRepositionVelocities", FLOAT, 1, UNIT_CONV_SPEED, 0.f, bparts,
last_repos_vel,
"Physical speeds at which the black holes repositioned most recently. "
"This is 0 for black holes that have never repositioned, or if the "
"simulation has been run without prescribed repositioning speed.");
list[41] = io_make_output_field(
"EddingtonFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
eddington_fraction,
"Accretion rates of black holes in units of their Eddington rates. "
"This is based on the unlimited accretion rates, so these fractions "
"can be above unity.");
list[42] = io_make_output_field(
"EnergyReservoirThresholds", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, bparts,
num_ngbs_to_heat,
"Minimum energy reservoir required for the black holes to do feedback, "
"expressed in units of the (constant) target heating temperature "
"increase.");
list[43] = io_make_output_field(
"BirthGasDensities", FLOAT, 1, UNIT_CONV_DENSITY, 0.f, bparts,
formation_gas_density,
"Physical densities of the converted part at the time of birth. "
......
......@@ -48,6 +48,9 @@ struct bpart {
/*! Black hole mass */
float mass;
/*! Black hole mass at the start of each step, prior to any nibbling */
float mass_at_start_of_step;
/* Particle cutoff radius. */
float h;
......@@ -115,6 +118,9 @@ struct bpart {
* radius (calculated as j_gas / h_BH, where j is specific ang. mom.) */
float circular_velocity_gas[3];
/*! Specific angular momentum of the gas around the black hole */
float spec_angular_momentum_gas[3];
/*! Multiplicative factor for accretion rates, from Rosas-Guevara et al.
* (2015) angular momentum based accretion disc model */
float f_visc;
......@@ -148,12 +154,48 @@ struct bpart {
* lower potential than all eligible neighbours) */
int number_of_reposition_attempts;
/* Velocity of most recent reposition jump */
float last_repos_vel;
/*! Total number of time steps in which the black hole was active. */
int number_of_time_steps;
/*! Total (physical) angular momentum accumulated by swallowing particles */
float swallowed_angular_momentum[3];
/*! Metallicity of ambient gas */
float gas_metal_mass_fraction;
/*! Accretion boost factor */
float accretion_boost_factor;
/*! Total (physical) angular momentum accumulated from subgrid accretion */
float accreted_angular_momentum[3];
/*! Cumulative target probability of heating any one neighbour particle */
float cumulative_target_prob;
/*! Cumulative actual probability of heating any one neighbour particle */
float cumulative_actual_prob;
/*! Instantaneous target probability of heating any one neighbour particle */
float target_heating_prob;
/*! Instantaneous coupling efficiency for feedback energy */
float epsilon_f;
/*! Cumulative coupling efficiency for feedback energy */
float cumulative_epsilon_f;
/*! Instantaneous heating temperature increase for feedback */
float AGN_delta_T;
/*! Instantaneous energy reservoir threshold (num-to-heat) */
float num_ngbs_to_heat;
/*! Eddington fractions */
float eddington_fraction;
/*! Union for the last high Eddington ratio point in time */
union {
......
......@@ -383,6 +383,52 @@ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
p_data->iron_mass_fraction_from_SNIa * gas_mass;
}
/**
* @brief Transfer chemistry data of a gas particle to a black hole.
*
* In contrast to `chemistry_add_part_to_bpart`, only a fraction of the
* masses stored in the gas particle are transferred here. Absolute masses
* of the gas particle are adjusted as well.
* Black holes don't store fractions so we need to add element masses.
*
* @param bp_data The black hole data to add to.
* @param p_data The gas data to use.
* @param nibble_mass The mass to be removed from the gas particle.
* @param nibble_fraction The fraction of the (original) mass of the gas
* particle that is removed.
*/
__attribute__((always_inline)) INLINE static void
chemistry_transfer_part_to_bpart(struct chemistry_bpart_data* bp_data,
struct chemistry_part_data* p_data,
const double nibble_mass,
const double nibble_fraction) {
bp_data->metal_mass_total += p_data->metal_mass_fraction_total * nibble_mass;
for (int i = 0; i < chemistry_element_count; ++i)
bp_data->metal_mass[i] += p_data->metal_mass_fraction[i] * nibble_mass;
bp_data->mass_from_SNIa += p_data->mass_from_SNIa * nibble_fraction;
bp_data->mass_from_SNII += p_data->mass_from_SNII * nibble_fraction;
bp_data->mass_from_AGB += p_data->mass_from_AGB * nibble_fraction;
/* Absolute masses, so need to reduce the gas particle */
p_data->mass_from_SNIa -= p_data->mass_from_SNIa * nibble_fraction;
p_data->mass_from_SNII -= p_data->mass_from_SNII * nibble_fraction;
p_data->mass_from_AGB -= p_data->mass_from_AGB * nibble_fraction;
bp_data->metal_mass_from_SNIa +=
p_data->metal_mass_fraction_from_SNIa * nibble_mass;
bp_data->metal_mass_from_SNII +=
p_data->metal_mass_fraction_from_SNII * nibble_mass;
bp_data->metal_mass_from_AGB +=
p_data->metal_mass_fraction_from_AGB * nibble_mass;
bp_data->iron_mass_from_SNIa +=
p_data->iron_mass_fraction_from_SNIa * nibble_mass;
/* TODO: Need to think about, and possibly implement, other transfers
* involving e.g. mass-weighted redshift of Fe enrichment */
}
/**
* @brief Add the chemistry data of a black hole to another one.
*
......@@ -510,4 +556,22 @@ chemistry_get_metal_mass_fraction_for_star_formation(
return p->chemistry_data.smoothed_metal_mass_fraction;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used for computing black hole ambient metallicities.
*
* Note that this is not smoothed, because we are already smoothing the
* individual particle metallicities within the BH kernel.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_metal_mass_fraction_for_black_holes(
const struct part* restrict p) {
return p->chemistry_data.metal_mass_fraction_total;
}
#endif /* SWIFT_CHEMISTRY_EAGLE_H */
......@@ -339,6 +339,30 @@ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
error("Loic: to be implemented");
}
/**
* @brief Transfer chemistry data of a gas particle to a black hole.
*
* In contrast to `chemistry_add_part_to_bpart`, only a fraction of the
* masses stored in the gas particle are transferred here. Absolute masses
* of the gas particle are adjusted as well.
* Black holes don't store fractions so we need to add element masses.
*
* Nothing to do here.
*
* @param bp_data The black hole data to add to.
* @param p_data The gas data to use.
* @param nibble_mass The mass to be removed from the gas particle.
* @param nibble_fraction The fraction of the (original) mass of the gas
* particle that is removed.
*/
__attribute__((always_inline)) INLINE static void
chemistry_transfer_part_to_bpart(struct chemistry_bpart_data* bp_data,
struct chemistry_part_data* p_data,
const double nibble_mass,
const double nibble_fraction) {
error("To be implemented.");
}
/**
* @brief Add the chemistry data of a black hole to another one.
*
......
......@@ -210,6 +210,23 @@ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
struct chemistry_bpart_data* bp_data,
const struct chemistry_part_data* p_data, const double gas_mass) {}
/**
* @brief Transfer chemistry data of a gas particle to a black hole.
*
* Nothing to do here.
*
* @param bp_data The black hole data to add to.
* @param p_data The gas data to use.
* @param nibble_mass The mass to be removed from the gas particle.
* @param nibble_fraction The fraction of the (original) mass of the gas
* particle that is removed.
*/
__attribute__((always_inline)) INLINE static void
chemistry_transfer_part_to_bpart(struct chemistry_bpart_data* bp_data,
struct chemistry_part_data* p_data,
const double nibble_mass,
const double nibble_fraction) {}
/**
* @brief Add the chemistry data of a black hole to another one.
*
......
......@@ -202,6 +202,23 @@ __attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
struct chemistry_bpart_data* bp_data,
const struct chemistry_part_data* p_data, const double gas_mass) {}
/**
* @brief Transfer chemistry data of a gas particle to a black hole.
*
* Nothing to do here.
*
* @param bp_data The black hole data to add to.
* @param p_data The gas data to use.
* @param nibble_mass The mass to be removed from the gas particle.
* @param nibble_fraction The fraction of the (original) mass of the gas
* particle that is removed.
*/
__attribute__((always_inline)) INLINE static void
chemistry_transfer_part_to_bpart(struct chemistry_bpart_data* bp_data,
struct chemistry_part_data* p_data,
const double nibble_mass,
const double nibble_fraction) {}
/**
* @brief Add the chemistry data of a black hole to another one.
*
......
This diff is collapsed.
......@@ -23,12 +23,16 @@
#include "error.h"
#include "feedback_properties.h"
#include "hydro_properties.h"
#include "star_formation.h"
#include "part.h"
#include "units.h"
#include <strings.h>
void compute_stellar_evolution(const struct feedback_props* feedback_props,
const struct star_formation* starform_props,
const struct hydro_props* hydro_props,
const struct phys_const* phys_const,
const struct cosmology* cosmo, struct spart* sp,
const struct unit_system* us, const double age,
const double dt);
......@@ -72,6 +76,8 @@ __attribute__((always_inline)) INLINE static void feedback_init_spart(
sp->feedback_data.to_collect.enrichment_weight_inv = 0.f;
sp->feedback_data.to_collect.ngb_mass = 0.f;
sp->feedback_data.to_collect.num_ngbs = 0;
sp->feedback_data.to_collect.gas_density = 0.f;
}
/**
......@@ -184,7 +190,8 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
const struct cosmology* cosmo, const struct unit_system* us,
const struct phys_const* phys_const, const double star_age_beg_step,
const double dt, const double time, const integertime_t ti_begin,
const int with_cosmology) {
const int with_cosmology, const struct star_formation* starform_props,
const struct hydro_props* hydro_props) {
#ifdef SWIFT_DEBUG_CHECKS
if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
......@@ -192,8 +199,9 @@ __attribute__((always_inline)) INLINE static void feedback_evolve_spart(
/* Compute amount of enrichment and feedback that needs to be done in this
* step */
compute_stellar_evolution(feedback_props, cosmo, sp, us, star_age_beg_step,
dt);
compute_stellar_evolution(feedback_props, starform_props, hydro_props,
phys_const, cosmo,
sp, us, star_age_beg_step, dt);
/* Decrease star mass by amount of mass distributed to gas neighbours */
sp->mass -= sp->feedback_data.to_distribute.mass;
......
......@@ -59,6 +59,8 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
/* Add mass of pj to neighbour mass of si */
si->feedback_data.to_collect.ngb_mass += mj;
si->feedback_data.to_collect.num_ngbs++;
si->feedback_data.to_collect.gas_density += mj * wi;
/* Add contribution of pj to normalisation of density weighted fraction
* which determines how much mass to distribute to neighbouring
......
......@@ -235,8 +235,30 @@ struct feedback_props {
/*! Wind delay time for SNII when using a fixed delay */
double SNII_wind_delay;
/*! Temperature increase induced by SNe feedback */
float SNe_deltaT_desired;
/*! Use variable temperature increase? */
int SNII_use_variable_delta_T;
/*! Buffer factor for numerical efficiency temperature */
double SNII_T_crit_factor;
/*! Should we use the instantaneous or birth density for determining dT? */
int SNII_use_instantaneous_density_for_dT;
/*! Factor to reduce sampling requirement if n_birth > n_birth_limit */
double SNII_sampling_nH_reduction_factor;
/*! Number of neighbours that should be heatable by SNII */
double SNII_delta_T_num_ngb_to_heat;
/*! Limiting number of neighbours that must be heatable by SNII */
double SNII_delta_T_num_ngb_to_heat_limit;
/*! Maximum temperature increase induced by SNII feedback [Kelvin] */
double SNII_delta_T_max;
double SNII_delta_T_min;
/*! (Constant) temperature increase induced by SNe feedback */
double SNe_deltaT_desired;
/*! Energy released by one supernova type II in cgs units */
double E_SNII_cgs;
......@@ -250,6 +272,15 @@ struct feedback_props {
/*! Maximal energy fraction for supernova type II feedback */
double f_E_max;
/*! Activate divergence-dependent SNII energy boost? */
int with_SNII_divergence_boost;
/*! Normalisation divergence for SNII energy boost */
double SNII_divergence_norm;
/*! Exponent for SNII energy dependence on divergence */
double SNII_divergence_exponent;
/*! Pivot point for the metallicity dependance of the feedback energy fraction
* model */
double Z_0;
......
......@@ -44,6 +44,12 @@ struct feedback_spart_data {
/*! Total mass (unweighted) of neighbouring gas particles */
float ngb_mass;
/*! Integer number of neighbouring gas particles */
int num_ngbs;
/*! Density of gas around the star particles */
float gas_density;
} to_collect;
/**
......
......@@ -431,6 +431,17 @@ __attribute__((always_inline)) INLINE static void hydro_set_viscosity_alpha(
p->viscosity.alpha = alpha;
}
/**
* @brief Get the value of the velocity divergence.
*
* @param p the particle of interest
* @param cosmo The cosmology model.
*/
__attribute__((always_inline)) INLINE static float
hydro_get_velocity_divergence(const struct part *restrict p) {
return p->viscosity.div_v;
}
/**
* @brief Update the value of the diffusive coefficients to the
* feedback reset value for the scheme.
......
......@@ -63,6 +63,8 @@ void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
struct part *parts = c->hydro.parts;
struct xpart *xparts = c->hydro.xparts;
const struct black_holes_props *props = e->black_holes_properties;
/* Early abort?
* (We only want cells for which we drifted the gas as these are
* the only ones that could have gas particles that have been flagged
......@@ -94,6 +96,13 @@ void runner_do_gas_swallow(struct runner *r, struct cell *c, int timer) {
/* Ignore inhibited particles (they have already been removed!) */
if (part_is_inhibited(p, e)) continue;
/* Update mass of associated gpart, to reflect potential changes from
* nibbling. In this case, we are already done. */
if (props->use_nibbling) {
p->gpart->mass = p->mass;
continue;
}
/* Get the ID of the black holes that will swallow this part */
const long long swallow_id =
black_holes_get_part_swallow_id(&p->black_holes_data);
......@@ -316,6 +325,10 @@ void runner_do_bh_swallow(struct runner *r, struct cell *c, int timer) {
/* Ignore inhibited particles (they have already been removed!) */
if (bpart_is_inhibited(cell_bp, e)) continue;
/* Update mass of associated gpart, to reflect potential changes from
* nibbling. */
if (props->use_nibbling) cell_bp->gpart->mass = cell_bp->mass;
/* Get the ID of the black holes that will swallow this bpart */
const long long swallow_id =
black_holes_get_bpart_swallow_id(&cell_bp->merger_data);
......
......@@ -75,9 +75,11 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct unit_system *us = e->internal_units;
const struct phys_const *phys_const = e->physical_constants;
const struct star_formation* starform_props = e->star_formation;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const struct cosmology *cosmo = e->cosmology;
const struct feedback_props *feedback_props = e->feedback_props;
const struct hydro_props* hydro_props = e->hydro_properties;
const float stars_h_max = e->hydro_properties->h_max;
const float stars_h_min = e->hydro_properties->h_min;
const float eps = e->stars_properties->h_tolerance;
......@@ -242,7 +244,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Compute the stellar evolution */
feedback_evolve_spart(sp, feedback_props, cosmo, us, phys_const,
star_age_beg_of_step, dt_enrichment,
e->time, ti_begin, with_cosmology);
e->time, ti_begin, with_cosmology,
starform_props, hydro_props);
} else {
/* Reset the feedback fields of the star particle */
......@@ -385,7 +388,8 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
/* Compute the stellar evolution */
feedback_evolve_spart(sp, feedback_props, cosmo, us, phys_const,
star_age_beg_of_step, dt_enrichment, e->time,
ti_begin, with_cosmology);
ti_begin, with_cosmology, starform_props,
hydro_props);
} else {
/* Reset the feedback fields of the star particle */
......
......@@ -464,10 +464,20 @@ INLINE static void star_formation_copy_properties(
sp->sf_data.birth_temperature = cooling_get_temperature(
phys_const, hydro_props, us, cosmo, cooling, p, xp);
sp->sf_data.birth_div_v = hydro_get_velocity_divergence(p);
sp->sf_data.birth_star_formation_rate = xp->sf_data.SFR;
/* Flag that this particle has not done feedback yet */
sp->f_E = -1.f;
sp->f_E_divergence_boost = -FLT_MAX;
sp->last_enrichment_time = sp->birth_time;
sp->count_since_last_enrichment = -1;
sp->delta_T_min = FLT_MAX;
sp->delta_T_max = -FLT_MAX;
sp->T_critical_fraction_min = FLT_MAX;
sp->T_critical_fraction_max = -FLT_MAX;
sp->T_sampling_fraction_min = FLT_MAX;
sp->T_sampling_fraction_max = -FLT_MAX;
}
/**
......
......@@ -73,6 +73,20 @@ star_formation_write_sparticles(const struct spart* sparts,
"Temperatures at the time of birth of the gas "
"particles that turned into stars");
return 2;
list[2] =
io_make_output_field("BirthVelocityDivergences", FLOAT, 1,
UNIT_CONV_FREQUENCY,
0.f, sparts, sf_data.birth_div_v,
"Velocity divergences at the time of birth of the "
"gas particles that turned into stars");
list[3] =
io_make_output_field(
"StarFormationRates", FLOAT, 1, UNIT_CONV_SFR, 0.f, sparts,
sf_data.birth_star_formation_rate,
"Star formation rates of the parent gas particle at the point where "
"they were converted to stars.");
return 4;
}
#endif /* SWIFT_STAR_FORMATION_EAGLE_IO_H */
......@@ -45,6 +45,13 @@ struct star_formation_spart_data {
/*! The birth temperature */
float birth_temperature;
/*! Velocity divergence of gas at birth */
float birth_div_v;
/*! Star formation rate of the parent gas particle */
float birth_star_formation_rate;
};
#endif /* SWIFT_EAGLE_STAR_FORMATION_STRUCT_H */
......@@ -68,8 +68,17 @@ __attribute__((always_inline)) INLINE static void stars_first_init_spart(
const int with_cosmology, const double scale_factor, const double time) {
sp->time_bin = 0;
sp->sf_data.birth_density = 0.f;
sp->sf_data.birth_density = -1.f;
sp->sf_data.birth_temperature = -1.f;
sp->sf_data.birth_div_v = FLT_MAX;
sp->f_E = -1.f;
sp->delta_T_min = FLT_MAX;
sp->delta_T_max = -FLT_MAX;
sp->T_critical_fraction_min = FLT_MAX;
sp->T_critical_fraction_max = -FLT_MAX;
sp->T_sampling_fraction_min = FLT_MAX;
sp->T_sampling_fraction_max = -FLT_MAX;
sp->f_E_divergence_boost = -FLT_MAX;
sp->count_since_last_enrichment = -1;
if (stars_properties->overwrite_birth_time)
......
......@@ -115,7 +115,7 @@ INLINE static void stars_write_particles(const struct spart *sparts,
const int with_cosmology) {
/* Say how much we want to write */
*num_fields = 8;
*num_fields = 15;
/* List what we want to write */
list[0] = io_make_output_field_convert_spart(
......@@ -158,6 +158,47 @@ INLINE static void stars_write_particles(const struct spart *sparts,
"FeedbackEnergyFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, sparts, f_E,
"Fractions of the canonical feedback energy that was used for the stars' "
"SNII feedback events");
list[8] = io_make_output_field(
"FeedbackMinimumDeltaT", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, sparts,
delta_T_min,
"Minimum temperature increase induced by SNII feedback from the stars.");
list[9] = io_make_output_field(
"FeedbackMaximumDeltaT", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, sparts,
delta_T_max,
"Maximum temperature increase induced by SNII feedback from the stars.");
list[10] = io_make_output_field(
"FeedbackMinimumCriticalFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
sparts, T_critical_fraction_min,
"Minimum temperature increase induced by SNII feedback from the stars in "
"units of the critical temperature for numerical efficiency.");
list[11] = io_make_output_field(
"FeedbackMaximumCriticalFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
sparts, T_critical_fraction_max,
"Maximum temperature increase induced by SNII feedback from the stars in "
"units of the critical temperature for numerical efficiency.");
list[12] = io_make_output_field(
"FeedbackMinimumSamplingFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
sparts, T_sampling_fraction_min,
"Minimum number of particles expected to be heated by each star "
"particle.");
list[13] = io_make_output_field(
"FeedbackMaximumSamplingFractions", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
sparts, T_sampling_fraction_max,
"Maximum number of particles expected to be heated by each star "
"particle.");
list[14] = io_make_output_field(
"FeedbackEnergyDivergenceBoosts", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f,
sparts, f_E_divergence_boost,
"Boost factor applied to energy feedback fractions because of negative "
"velocity divergence at birth of the star particles. This is already "
"incorporated in the value of FeedbackEnergyFractions.");
}
/**
......
......@@ -89,6 +89,21 @@ struct spart {
/*! Feedback energy fraction */
float f_E;
/*! Value of divergence boost factor (within f_E) */
float f_E_divergence_boost;
/*! Feedback temperature increase */
float delta_T_min;
float delta_T_max;
/*! Heating temperature in units of critical temperature */
float T_critical_fraction_min;
float T_critical_fraction_max;
/*! Heating temperature in units of T to heat one particle */
float T_sampling_fraction_min;
float T_sampling_fraction_max;
/*! Star formation struct */
struct star_formation_spart_data sf_data;
......
......@@ -63,6 +63,10 @@ tot_num_parts = tot_num_parts.astype(np.int64)
for i in range(6):
tot_num_parts[i] += (np.int64(tot_num_parts_high_word[i]) << 32)
tot_num_parts_swift = np.copy(tot_num_parts)
tot_num_parts_swift[2] += tot_num_parts_swift[3]
tot_num_parts_swift[3] = 0
# Some basic information
print("Reading", tot_num_parts, "particles from", num_files, "files.")
......@@ -85,9 +89,9 @@ output_file = h5.File(output_file_name, "w-")
# Header
grp = output_file.create_group("/Header")
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["NumPart_Total"] = tot_num_parts
grp.attrs["NumPart_Total"] = tot_num_parts_swift
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = tot_num_parts
grp.attrs["NumPart_ThisFile"] = tot_num_parts_swift
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["BoxSize"] = box_size
grp.attrs["Flag_Entropy_ICs"] = entropy_flag
......@@ -98,6 +102,8 @@ if tot_num_parts[0] > 0:
grp0 = output_file.create_group("/PartType0")
if tot_num_parts[1] > 0:
grp1 = output_file.create_group("/PartType1")
if tot_num_parts_swift[2] > 0:
grp2 = output_file.create_group("/PartType2")
if tot_num_parts[4] > 0:
grp4 = output_file.create_group("/PartType4")
if tot_num_parts[5] > 0:
......@@ -147,6 +153,12 @@ if tot_num_parts[1] > 0:
create_set(grp1, "Masses", tot_num_parts[1], 1, "f")
create_set(grp1, "ParticleIDs", tot_num_parts[1], 1, "l")
if tot_num_parts_swift[2] > 0:
create_set(grp2, "Coordinates", tot_num_parts_swift[2], 3, "d")
create_set(grp2, "Velocities", tot_num_parts_swift[2], 3, "f")
create_set(grp2, "Masses", tot_num_parts_swift[2], 1, "f")
create_set(grp2, "ParticleIDs", tot_num_parts_swift[2], 1, "l")
if tot_num_parts[4] > 0:
create_set(grp4, "Coordinates", tot_num_parts[4], 3, "d")
create_set(grp4, "Velocities", tot_num_parts[4], 3, "f")
......@@ -188,6 +200,8 @@ for f in range(num_files):
": num_parts = [",
num_parts[0],
num_parts[1],
num_parts[2],
num_parts[3],
num_parts[4],
num_parts[5],
"]",
......@@ -205,6 +219,14 @@ for f in range(num_files):
def copy_grp_same_name(name, ptype):
copy_grp(name, name, ptype)
def copy_grp_pt3(name):
full_name_new = "/PartType2/" + name
full_name_old = "/PartType3/" + name
output_file[full_name_new][
cumul_parts[2] : cumul_parts[2] + num_parts[3]
] = file[full_name_old]
if num_parts[0] > 0:
copy_grp_same_name("Coordinates", 0)
copy_grp_same_name("Velocities", 0)
......@@ -220,6 +242,21 @@ for f in range(num_files):
if DM_mass == 0.0: # Do not overwrite values if there was a mass table
copy_grp_same_name("Masses", 1)
if num_parts[2] > 0:
copy_grp_same_name("Coordinates", 2)
copy_grp_same_name("Velocities", 2)
copy_grp_same_name("ParticleIDs", 2)
copy_grp_same_name("Masses", 2)
# Need to update part counter for pt2 already here, so we append 3 in correct place
cumul_parts[2] += num_parts[2]
if num_parts[3] > 0:
copy_grp_pt3("Coordinates")
copy_grp_pt3("Velocities")
copy_grp_pt3("ParticleIDs")
copy_grp_pt3("Masses")
if num_parts[4] > 0:
copy_grp_same_name("Coordinates", 4)
copy_grp_same_name("Velocities", 4)
......@@ -234,6 +271,7 @@ for f in range(num_files):
cumul_parts[0] += num_parts[0]
cumul_parts[1] += num_parts[1]
cumul_parts[2] += num_parts[3] # Need to adjust for added pt-3
cumul_parts[4] += num_parts[4]
cumul_parts[5] += num_parts[5]
file.close()
......