Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • dc-oman1/swiftsim
  • swift/swiftsim
  • pdraper/swiftsim
  • tkchan/swiftsim
  • dc-turn5/swiftsim
5 results
Select Git revision
Show changes
...@@ -25,6 +25,8 @@ ...@@ -25,6 +25,8 @@
float initial_mass_function_get_exponent( float initial_mass_function_get_exponent(
const struct initial_mass_function *imf, float mass_min, float mass_max); const struct initial_mass_function *imf, float mass_min, float mass_max);
void initial_mass_function_print(const struct initial_mass_function *imf); void initial_mass_function_print(const struct initial_mass_function *imf);
float initial_mass_function_sample(const struct initial_mass_function *imf,
float f);
void initial_mass_function_integrate(const struct initial_mass_function *imf, void initial_mass_function_integrate(const struct initial_mass_function *imf,
float *data, size_t count, float *data, size_t count,
...@@ -35,8 +37,11 @@ float initial_mass_function_get_integral_xi( ...@@ -35,8 +37,11 @@ float initial_mass_function_get_integral_xi(
const struct initial_mass_function *imf, float m1, float m2); const struct initial_mass_function *imf, float m1, float m2);
float initial_mass_function_get_imf(const struct initial_mass_function *imf, float initial_mass_function_get_imf(const struct initial_mass_function *imf,
float m); float m);
float initial_mass_function_get_integral_imf( float initial_mass_function_get_imf_mass_fraction(
const struct initial_mass_function *imf, const float m1, const float m2); const struct initial_mass_function *imf, const float m1, const float m2);
float initial_mass_function_get_imf_number_fraction(
const struct initial_mass_function *imf, const float m1, const float m2);
void initial_mass_function_compute_coefficients( void initial_mass_function_compute_coefficients(
struct initial_mass_function *imf); struct initial_mass_function *imf);
...@@ -58,4 +63,12 @@ void initial_mass_function_restore(struct initial_mass_function *imf, ...@@ -58,4 +63,12 @@ void initial_mass_function_restore(struct initial_mass_function *imf,
const struct stellar_model *sm); const struct stellar_model *sm);
void initial_mass_function_clean(struct initial_mass_function *imf); void initial_mass_function_clean(struct initial_mass_function *imf);
double initial_mass_function_sample_power_law(double min_mass, double max_mass,
double exp, double x);
void initial_mass_function_compute_Mc_Md_Mtot(
const struct initial_mass_function *imf, double *M_continuous,
double *M_discrete, double *M_tot);
#endif // SWIFT_INITIAL_MASS_FUNCTION_GEAR_H #endif // SWIFT_INITIAL_MASS_FUNCTION_GEAR_H
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#include "stellar_evolution.h" #include "stellar_evolution.h"
/* Include local headers */ /* Include local headers */
#include "exp10.h"
#include "hdf5_functions.h" #include "hdf5_functions.h"
#include "initial_mass_function.h" #include "initial_mass_function.h"
#include "lifetime.h" #include "lifetime.h"
...@@ -32,6 +33,8 @@ ...@@ -32,6 +33,8 @@
#include <math.h> #include <math.h>
#include <stddef.h> #include <stddef.h>
#define DEFAULT_STAR_MINIMAL_GRAVITY_MASS_MSUN 1e-1
/** /**
* @brief Print the stellar model. * @brief Print the stellar model.
* *
...@@ -80,6 +83,81 @@ int stellar_evolution_compute_integer_number_supernovae( ...@@ -80,6 +83,81 @@ int stellar_evolution_compute_integer_number_supernovae(
return number_supernovae_i + (rand_sn < frac_sn); return number_supernovae_i + (rand_sn < frac_sn);
} }
/**
* @brief Update the #spart mass from the supernovae ejected mass.
*
* This function deals with each star_type.
*
* Note: This function is called by
* stellar_evolution_compute_discrete_feedback_properties() and
* stellar_evolution_compute_continuous_feedback_properties().
*
* @param sp The particle to act upon
* @param sm The #stellar_model structure.
*/
void stellar_evolution_sn_apply_ejected_mass(struct spart* restrict sp,
const struct stellar_model* sm) {
/* If a star is a discrete star */
if (sp->star_type == single_star) {
const int null_mass = (sp->mass == sp->feedback_data.mass_ejected);
const int negative_mass = (sp->mass < sp->feedback_data.mass_ejected);
if (null_mass) {
message("Star %lld (m_star = %e, m_ej = %e) completely exploded!", sp->id,
sp->mass, sp->feedback_data.mass_ejected);
/* If the star ejects all its mass (for very massive stars), give it a
zero mass so that we know it has exploded.
We do not remove the star from the simulation to keep track of its
properties, e.g. to check the IMF sampling (with sinks).
Bug fix (28.04.2024): The mass of the star should not be set to
0 because of gravity. So, we give some minimal value. */
sp->mass = sm->discrete_star_minimal_gravity_mass;
/* If somehow the star has a negative mass, we have a problem. */
} else if (negative_mass) {
error(
"(Discrete star) Negative mass (m_star = %e, m_ej = %e), skipping "
"current star: %lli",
sp->mass, sp->feedback_data.mass_ejected, sp->id);
/* Reset everything */
sp->feedback_data.number_snia = 0;
sp->feedback_data.number_snii = 0;
sp->feedback_data.mass_ejected = 0;
/* Reset energy to avoid injecting anything in the
runner_iact_nonsym_feedback_apply() */
sp->feedback_data.energy_ejected = 0;
return;
} else {
/* Update the mass */
sp->mass -= sp->feedback_data.mass_ejected;
}
/* If the star is the continuous part of the IMF or the enteire IMF */
} else {
/* Check if we can ejected the required amount of elements. */
const int negative_mass = (sp->mass <= sp->feedback_data.mass_ejected);
if (negative_mass) {
warning(
"(Continuous star) Negative mass (m_star = %e, m_ej = %e), skipping "
"current star: %lli",
sp->mass, sp->feedback_data.mass_ejected, sp->id);
/* Reset everything */
sp->feedback_data.number_snia = 0;
sp->feedback_data.number_snii = 0;
sp->feedback_data.mass_ejected = 0;
/* Reset energy to avoid injecting anything in the
runner_iact_nonsym_feedback_apply() */
sp->feedback_data.energy_ejected = 0;
return;
}
/* Update the mass */
sp->mass -= sp->feedback_data.mass_ejected;
}
}
/** /**
* @brief Compute the feedback properties. * @brief Compute the feedback properties.
* *
...@@ -121,19 +199,8 @@ void stellar_evolution_compute_continuous_feedback_properties( ...@@ -121,19 +199,8 @@ void stellar_evolution_compute_continuous_feedback_properties(
sp->feedback_data.mass_ejected = mass_frac_snii * sp->sf_data.birth_mass + sp->feedback_data.mass_ejected = mass_frac_snii * sp->sf_data.birth_mass +
mass_snia * phys_const->const_solar_mass; mass_snia * phys_const->const_solar_mass;
/* Check if we can ejected the required amount of elements. */ /* Removes the ejected mass from the star */
const int negative_mass = sp->mass <= sp->feedback_data.mass_ejected; stellar_evolution_sn_apply_ejected_mass(sp, sm);
if (negative_mass) {
message("Negative mass, skipping current star: %lli", sp->id);
/* Reset everything */
sp->feedback_data.number_snia = 0;
sp->feedback_data.number_snii = 0;
sp->feedback_data.mass_ejected = 0;
return;
}
/* Update the mass */
sp->mass -= sp->feedback_data.mass_ejected;
/* Now deal with the metals */ /* Now deal with the metals */
...@@ -214,19 +281,8 @@ void stellar_evolution_compute_discrete_feedback_properties( ...@@ -214,19 +281,8 @@ void stellar_evolution_compute_discrete_feedback_properties(
/* Transform into internal units */ /* Transform into internal units */
sp->feedback_data.mass_ejected *= phys_const->const_solar_mass; sp->feedback_data.mass_ejected *= phys_const->const_solar_mass;
/* Check if we can ejected the required amount of elements. */ /* Removes the ejected mass from the star */
const int negative_mass = sp->mass <= sp->feedback_data.mass_ejected; stellar_evolution_sn_apply_ejected_mass(sp, sm);
if (negative_mass) {
message("Negative mass, skipping current star: %lli", sp->id);
/* Reset everything */
sp->feedback_data.number_snia = 0;
sp->feedback_data.number_snii = 0;
sp->feedback_data.mass_ejected = 0;
return;
}
/* Update the mass */
sp->mass -= sp->feedback_data.mass_ejected;
/* Get the SNIa yields */ /* Get the SNIa yields */
const float* snia_yields = supernovae_ia_get_yields(&sm->snia); const float* snia_yields = supernovae_ia_get_yields(&sm->snia);
...@@ -262,6 +318,110 @@ void stellar_evolution_compute_discrete_feedback_properties( ...@@ -262,6 +318,110 @@ void stellar_evolution_compute_discrete_feedback_properties(
} }
} }
/**
* @brief Evolve an individual star represented by a #spart.
*
* This function compute the SN rate and yields before sending
* this information to a different MPI rank.
* It also compute the supernovae energy to be released by the
* star.
*
* Here I am using Myr-solar mass units internally in order to
* avoid numerical errors.
*
* Note: This function treats the case of single/individual stars.
*
* @param sp The particle to act upon
* @param sm The #stellar_model structure.
* @param cosmo The current cosmological model.
* @param us The unit system.
* @param phys_const The physical constants in the internal unit system.
* @param ti_begin The #integertime_t at the begining of the step.
* @param star_age_beg_step The age of the star at the star of the time-step in
* internal units.
* @param dt The time-step size of this star in internal units.
*/
void stellar_evolution_evolve_individual_star(
struct spart* restrict sp, const struct stellar_model* sm,
const struct cosmology* cosmo, const struct unit_system* us,
const struct phys_const* phys_const, const integertime_t ti_begin,
const double star_age_beg_step, const double dt) {
/* Check that this function is called for single_star only. */
if (sp->star_type != single_star) {
error("This function can only be called for single/individual star!");
}
/* Convert the inputs */
const double conversion_to_myr = phys_const->const_year * 1e6;
const double star_age_end_step_myr =
(star_age_beg_step + dt) / conversion_to_myr;
const double star_age_beg_step_myr = star_age_beg_step / conversion_to_myr;
/* Get the metallicity */
const float metallicity =
chemistry_get_star_total_metal_mass_fraction_for_feedback(sp);
const float log_mass = log10(sp->mass / phys_const->const_solar_mass);
const float lifetime_myr = pow(10, lifetime_get_log_lifetime_from_mass(
&sm->lifetime, log_mass, metallicity));
/* if the lifetime is outside the interval */
if ((lifetime_myr < star_age_beg_step_myr) ||
(lifetime_myr > star_age_end_step_myr))
return;
message(
"(%lld) lifetime_myr=%g %g star_age_beg_step=%g star_age_end_step=%g "
"(%g)",
sp->id, lifetime_myr, lifetime_myr * conversion_to_myr,
star_age_beg_step_myr, star_age_end_step_myr,
sp->mass / phys_const->const_solar_mass);
/* This is needed by stellar_evolution_compute_discrete_feedback_properties(),
but this is not used inside the function. */
const float m_init = 0;
/* Get the integer number of supernovae */
const int number_snia = 0;
const int number_snii = 1;
/* Save the number of supernovae */
sp->feedback_data.number_snia = 0;
sp->feedback_data.number_snii = number_snii;
/* this is needed for stellar_evolution_compute_discrete_feedback_properties
*/
const float m_beg_step = sp->mass / phys_const->const_solar_mass;
const float m_end_step = sp->mass / phys_const->const_solar_mass;
const float m_avg = 0.5 * (m_beg_step + m_end_step);
/* Compute the yields */
stellar_evolution_compute_discrete_feedback_properties(
sp, sm, phys_const, m_beg_step, m_end_step, m_init, number_snia,
number_snii);
/* Compute the supernovae energy associated to the stellar particle */
const float energy_conversion =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) / 1e51;
/* initialize */
sp->feedback_data.energy_ejected = 0;
/* snia contribution */
const float snia_energy = sm->snia.energy_per_supernovae;
sp->feedback_data.energy_ejected +=
sp->feedback_data.number_snia * snia_energy;
/* snii contribution */
const float snii_energy =
supernovae_ii_get_energy_from_progenitor_mass(&sm->snii, m_avg) /
energy_conversion;
sp->feedback_data.energy_ejected +=
sp->feedback_data.number_snii * snii_energy;
}
/** /**
* @brief Evolve the stellar properties of a #spart. * @brief Evolve the stellar properties of a #spart.
* *
...@@ -273,6 +433,10 @@ void stellar_evolution_compute_discrete_feedback_properties( ...@@ -273,6 +433,10 @@ void stellar_evolution_compute_discrete_feedback_properties(
* Here I am using Myr-solar mass units internally in order to * Here I am using Myr-solar mass units internally in order to
* avoid numerical errors. * avoid numerical errors.
* *
* Note: This function treats the case of particles representing the whole IMF
* (star_type = star_population) and the particles representing only the
* continuous part of the IMF (star_type = star_population_continuous_IMF).
*
* @param sp The particle to act upon * @param sp The particle to act upon
* @param sm The #stellar_model structure. * @param sm The #stellar_model structure.
* @param cosmo The current cosmological model. * @param cosmo The current cosmological model.
...@@ -289,6 +453,14 @@ void stellar_evolution_evolve_spart( ...@@ -289,6 +453,14 @@ void stellar_evolution_evolve_spart(
const struct phys_const* phys_const, const integertime_t ti_begin, const struct phys_const* phys_const, const integertime_t ti_begin,
const double star_age_beg_step, const double dt) { const double star_age_beg_step, const double dt) {
/* Check that this function is called for populations of stars and not
individual stars. */
if (sp->star_type == single_star) {
error(
"This function can only be called for sparts representing stars "
"populations!");
}
/* Convert the inputs */ /* Convert the inputs */
const double conversion_to_myr = phys_const->const_year * 1e6; const double conversion_to_myr = phys_const->const_year * 1e6;
const double star_age_beg_step_myr = star_age_beg_step / conversion_to_myr; const double star_age_beg_step_myr = star_age_beg_step / conversion_to_myr;
...@@ -318,6 +490,30 @@ void stellar_evolution_evolve_spart( ...@@ -318,6 +490,30 @@ void stellar_evolution_evolve_spart(
*/ */
if (m_end_step >= m_beg_step) return; if (m_end_step >= m_beg_step) return;
/* Star particles representing only the continuous part of the IMF need a
special treatment. They do not contain stars above the mass that separate the
IMF into two parts (variable called minimal_discrete_mass_Msun in the sink
module). So, if m_end_step > minimal_discrete_mass_Msun, you don't do
feedback. Note that the sm structure contains different information for the
'first stars' and the 'late stars'. The right sm data is passed to this
function so we do not need any special treatment here. */
if (sp->star_type == star_population_continuous_IMF) {
/* If it's not time yet for feedback, exit. Notice that both masses are in
solar mass. */
if (m_end_step > sm->imf.minimal_discrete_mass_Msun) {
return;
}
/* If we are in a case where
m_beg_step > minimal_discrete_mass_Msun > m_end_step,
then we need to be careful. We don't want feedback from the discrete
part, only the continuous part. Hence, we need to update m_beg_step.
*/
if (m_beg_step > sm->imf.minimal_discrete_mass_Msun) {
m_beg_step = sm->imf.minimal_discrete_mass_Msun;
}
}
/* Check if the star can produce a supernovae */ /* Check if the star can produce a supernovae */
const int can_produce_snia = const int can_produce_snia =
supernovae_ia_can_explode(&sm->snia, m_end_step, m_beg_step); supernovae_ia_can_explode(&sm->snia, m_end_step, m_beg_step);
...@@ -327,8 +523,14 @@ void stellar_evolution_evolve_spart( ...@@ -327,8 +523,14 @@ void stellar_evolution_evolve_spart(
/* Is it possible to generate a supernovae? */ /* Is it possible to generate a supernovae? */
if (!can_produce_snia && !can_produce_snii) return; if (!can_produce_snia && !can_produce_snii) return;
/* Compute the initial mass */ /* Compute the initial mass. The initial mass is different if the star
const float m_init = sp->sf_data.birth_mass / phys_const->const_solar_mass; particle is of type 'star_population' or
'star_population_continuous_IMF'. The function call treats both cases. */
const float m_init =
stellar_evolution_compute_initial_mass(sp, sm, phys_const);
/* Then, for 'star_population_continuous_IMF', everything remain the same as
with the "old" 'star_population'! */
/* Compute number of SNIa */ /* Compute number of SNIa */
float number_snia_f = 0; float number_snia_f = 0;
...@@ -434,6 +636,21 @@ int stellar_evolution_get_element_index(const struct stellar_model* sm, ...@@ -434,6 +636,21 @@ int stellar_evolution_get_element_index(const struct stellar_model* sm,
return -1; return -1;
} }
/**
* @brief Get the solar abundance of the element .
*
* @param sm The #stellar_model.
* @param element_name The element name.
*/
float stellar_evolution_get_solar_abundance(const struct stellar_model* sm,
const char* element_name) {
int element_index = stellar_evolution_get_element_index(sm, element_name);
float solar_abundance = sm->solar_abundances[element_index];
return solar_abundance;
}
/** /**
* @brief Read the name of all the elements present in the tables. * @brief Read the name of all the elements present in the tables.
* *
...@@ -487,6 +704,46 @@ void stellar_evolution_read_elements(struct stellar_model* sm, ...@@ -487,6 +704,46 @@ void stellar_evolution_read_elements(struct stellar_model* sm,
} }
} }
/**
* @brief Read the solar abundances.
*
* @param parameter_file The parsed parameter file.
* @param data The properties to initialise.
*/
void stellar_evolution_read_solar_abundances(struct stellar_model* sm,
struct swift_params* params) {
#if defined(HAVE_HDF5)
/* Get the yields table */
char filename[DESCRIPTION_BUFFER_SIZE];
parser_get_param_string(params, "GEARFeedback:yields_table", filename);
/* Open file. */
hid_t file_id = H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT);
if (file_id < 0) error("unable to open file %s.\n", filename);
/* Open group. */
hid_t group_id = H5Gopen(file_id, "Data", H5P_DEFAULT);
if (group_id < 0) error("unable to open group Data.\n");
/* Read the data */
io_read_array_attribute(group_id, "SolarMassAbundances", FLOAT,
sm->solar_abundances, GEAR_CHEMISTRY_ELEMENT_COUNT);
/* Close group */
hid_t status = H5Gclose(group_id);
if (status < 0) error("error closing group.");
/* Close file */
status = H5Fclose(file_id);
if (status < 0) error("error closing file.");
#else
message("Cannot read the solar abundances without HDF5");
#endif
}
/** /**
* @brief Initialize the global properties of the stellar evolution scheme. * @brief Initialize the global properties of the stellar evolution scheme.
* *
...@@ -505,6 +762,9 @@ void stellar_evolution_props_init(struct stellar_model* sm, ...@@ -505,6 +762,9 @@ void stellar_evolution_props_init(struct stellar_model* sm,
/* Read the list of elements */ /* Read the list of elements */
stellar_evolution_read_elements(sm, params); stellar_evolution_read_elements(sm, params);
/* Read the solar abundances */
stellar_evolution_read_solar_abundances(sm, params);
/* Use the discrete yields approach? */ /* Use the discrete yields approach? */
sm->discrete_yields = sm->discrete_yields =
parser_get_param_int(params, "GEARFeedback:discrete_yields"); parser_get_param_int(params, "GEARFeedback:discrete_yields");
...@@ -521,6 +781,20 @@ void stellar_evolution_props_init(struct stellar_model* sm, ...@@ -521,6 +781,20 @@ void stellar_evolution_props_init(struct stellar_model* sm,
/* Initialize the supernovae II model */ /* Initialize the supernovae II model */
supernovae_ii_init(&sm->snii, params, sm, us); supernovae_ii_init(&sm->snii, params, sm, us);
/* Initialize the minimal gravity mass for the stars */
/* const float default_star_minimal_gravity_mass_Msun = 1e-1; */
sm->discrete_star_minimal_gravity_mass = parser_get_opt_param_float(
params, "GEARFeedback:discrete_star_minimal_gravity_mass_Msun",
DEFAULT_STAR_MINIMAL_GRAVITY_MASS_MSUN);
/* Convert from M_sun to internal units */
sm->discrete_star_minimal_gravity_mass *= phys_const->const_solar_mass;
if (engine_rank == 0) {
message("discrete_star_minimal_gravity_mass: (internal units) %e",
sm->discrete_star_minimal_gravity_mass);
}
} }
/** /**
...@@ -585,3 +859,38 @@ void stellar_evolution_clean(struct stellar_model* sm) { ...@@ -585,3 +859,38 @@ void stellar_evolution_clean(struct stellar_model* sm) {
supernovae_ia_clean(&sm->snia); supernovae_ia_clean(&sm->snia);
supernovae_ii_clean(&sm->snii); supernovae_ii_clean(&sm->snii);
} }
/**
* @brief Computes the initial mass of a #spart. This function distinguishes
* between the stellar particle representing a whole IMF and the stellar
* particles representing only the continuous part.
*
* @param sp The particle for which we compute the initial mass.
* @param sm The #stellar_model structure.
* @param phys_const the physical constants in internal units.
* @param (return) m_init Initial mass of the star particle (in M_sun).
*/
float stellar_evolution_compute_initial_mass(
const struct spart* restrict sp, const struct stellar_model* sm,
const struct phys_const* phys_const) {
const struct initial_mass_function* imf = &sm->imf;
switch (sp->star_type) {
case star_population:
return sp->sf_data.birth_mass / phys_const->const_solar_mass;
case star_population_continuous_IMF: {
double M_IMF_tot, M_d_dummy, M_c_dummy;
initial_mass_function_compute_Mc_Md_Mtot(imf, &M_c_dummy, &M_d_dummy,
&M_IMF_tot);
/* No need to convert from internal units to M_sun because the masses are
already in solar masses (to avoid numerical errors) */
return M_IMF_tot;
}
case single_star:
return sp->sf_data.birth_mass / phys_const->const_solar_mass;
default: {
error("This star_type (%d) is not implemented!", sp->star_type);
return -1.0;
}
}
}
...@@ -46,6 +46,12 @@ void stellar_evolution_compute_discrete_feedback_properties( ...@@ -46,6 +46,12 @@ void stellar_evolution_compute_discrete_feedback_properties(
const float m_end_step, const float m_init, const int number_snia, const float m_end_step, const float m_init, const int number_snia,
const int number_snii); const int number_snii);
void stellar_evolution_evolve_individual_star(
struct spart* restrict sp, const struct stellar_model* sm,
const struct cosmology* cosmo, const struct unit_system* us,
const struct phys_const* phys_const, const integertime_t ti_begin,
const double star_age_beg_step, const double dt);
void stellar_evolution_evolve_spart( void stellar_evolution_evolve_spart(
struct spart* restrict sp, const struct stellar_model* sm, struct spart* restrict sp, const struct stellar_model* sm,
const struct cosmology* cosmo, const struct unit_system* us, const struct cosmology* cosmo, const struct unit_system* us,
...@@ -56,9 +62,12 @@ const char* stellar_evolution_get_element_name(const struct stellar_model* sm, ...@@ -56,9 +62,12 @@ const char* stellar_evolution_get_element_name(const struct stellar_model* sm,
int i); int i);
int stellar_evolution_get_element_index(const struct stellar_model* sm, int stellar_evolution_get_element_index(const struct stellar_model* sm,
const char* element_name); const char* element_name);
float stellar_evolution_get_solar_abundance(const struct stellar_model* sm,
const char* element_name);
void stellar_evolution_read_elements(struct stellar_model* sm, void stellar_evolution_read_elements(struct stellar_model* sm,
struct swift_params* params); struct swift_params* params);
void stellar_evolution_read_solar_abundances(struct stellar_model* sm,
struct swift_params* params);
void stellar_evolution_props_init(struct stellar_model* sm, void stellar_evolution_props_init(struct stellar_model* sm,
const struct phys_const* phys_const, const struct phys_const* phys_const,
const struct unit_system* us, const struct unit_system* us,
...@@ -70,4 +79,7 @@ void stellar_evolution_restore(struct stellar_model* sm, FILE* stream); ...@@ -70,4 +79,7 @@ void stellar_evolution_restore(struct stellar_model* sm, FILE* stream);
void stellar_evolution_clean(struct stellar_model* sm); void stellar_evolution_clean(struct stellar_model* sm);
float stellar_evolution_compute_initial_mass(
const struct spart* restrict sp, const struct stellar_model* sm,
const struct phys_const* phys_consts);
#endif // SWIFT_STELLAR_EVOLUTION_GEAR_H #endif // SWIFT_STELLAR_EVOLUTION_GEAR_H
...@@ -39,13 +39,17 @@ struct initial_mass_function { ...@@ -39,13 +39,17 @@ struct initial_mass_function {
/*! Mass limits between IMF parts (n_parts + 1 elements). */ /*! Mass limits between IMF parts (n_parts + 1 elements). */
float *mass_limits; float *mass_limits;
/*! Mass fraction computed at the interface between two IMF parts (n_parts + 1
* elements). */
float *mass_fraction;
/*! Exponent of each IMF parts (n_parts elements). */ /*! Exponent of each IMF parts (n_parts elements). */
float *exp; float *exp;
/*! Coefficient of each IMF parts (n_parts elements). */ /*! Coefficient of each IMF parts (n_parts elements). */
float *coef; float *coef;
/*! Number of parts in the function. */ /*! Number of parts (segments) in the function. */
int n_parts; int n_parts;
/*! Minimal mass contained in mass_limits, copied for more clarity. */ /*! Minimal mass contained in mass_limits, copied for more clarity. */
...@@ -53,6 +57,19 @@ struct initial_mass_function { ...@@ -53,6 +57,19 @@ struct initial_mass_function {
/*! Maximal mass contained in mass_limits, copied for more clarity. */ /*! Maximal mass contained in mass_limits, copied for more clarity. */
float mass_max; float mass_max;
/*! Total number of stars (per mass unit) in the IMF. */
float N_tot;
/*! Probability to generate a star out of the continuous part of the IMF. */
float sink_Pc;
/*! Stellar mass of the continous part of the IMF (in solar mass). */
float stellar_particle_mass_Msun;
/*! Minimal mass of stars represented by discrete particles (in solar mass).
*/
float minimal_discrete_mass_Msun;
}; };
/** /**
...@@ -174,6 +191,9 @@ struct stellar_model { ...@@ -174,6 +191,9 @@ struct stellar_model {
/*! Name of the different elements */ /*! Name of the different elements */
char elements_name[GEAR_CHEMISTRY_ELEMENT_COUNT * GEAR_LABELS_SIZE]; char elements_name[GEAR_CHEMISTRY_ELEMENT_COUNT * GEAR_LABELS_SIZE];
/* Solar mass abundances read from the chemistry table */
float solar_abundances[GEAR_CHEMISTRY_ELEMENT_COUNT];
/*! The initial mass function */ /*! The initial mass function */
struct initial_mass_function imf; struct initial_mass_function imf;
...@@ -191,6 +211,19 @@ struct stellar_model { ...@@ -191,6 +211,19 @@ struct stellar_model {
/* Filename of the yields table */ /* Filename of the yields table */
char yields_table[FILENAME_BUFFER_SIZE]; char yields_table[FILENAME_BUFFER_SIZE];
/* Minimal gravity mass after a discrete star has completely exploded.
This will be the mass of the gpart's friend of the star. The mass of the
star will be 0 after it losses all its mass.
The purpose of this is to avoid zero mass for the gravitsy
computations. We keep the star so that we know it *existed* and we can
extract its properties at the end of a run. If we remove the star, then
we do not have any information about its existence.
However, since the star is dead/inexistent, the gravity mass must be small
so that it does not drastically alter the dynamics of the systems. */
float discrete_star_minimal_gravity_mass;
}; };
#endif // SWIFT_STELLAR_EVOLUTION_STRUCT_GEAR_H #endif // SWIFT_STELLAR_EVOLUTION_STRUCT_GEAR_H
...@@ -192,7 +192,7 @@ void supernovae_ia_read_yields(struct supernovae_ia *snia, ...@@ -192,7 +192,7 @@ void supernovae_ia_read_yields(struct supernovae_ia *snia,
io_read_array_attribute(group_id, "data", FLOAT, yields, nval); io_read_array_attribute(group_id, "data", FLOAT, yields, nval);
/* Read the labels */ /* Read the labels */
char *labels = malloc(nval * GEAR_LABELS_SIZE); char *labels = (char *)calloc(nval, GEAR_LABELS_SIZE);
io_read_string_array_attribute(group_id, "elts", labels, nval, io_read_string_array_attribute(group_id, "elts", labels, nval,
GEAR_LABELS_SIZE); GEAR_LABELS_SIZE);
......
...@@ -65,7 +65,7 @@ int current_fof_attach_type; ...@@ -65,7 +65,7 @@ int current_fof_attach_type;
int current_fof_ignore_type; int current_fof_ignore_type;
/* Are we timing calculating group properties in the FOF? */ /* Are we timing calculating group properties in the FOF? */
//#define WITHOUT_GROUP_PROPS // #define WITHOUT_GROUP_PROPS
/** /**
* @brief Properties of a group used for black hole seeding * @brief Properties of a group used for black hole seeding
...@@ -709,6 +709,33 @@ __attribute__((always_inline)) INLINE static int gpart_is_ignorable( ...@@ -709,6 +709,33 @@ __attribute__((always_inline)) INLINE static int gpart_is_ignorable(
return current_fof_ignore_type & (1 << (gp->type + 1)); return current_fof_ignore_type & (1 << (gp->type + 1));
} }
/**
* @brief Returns whether a foreign #gpart is of the 'attachable' kind.
*/
__attribute__((always_inline)) INLINE static int gpart_foreign_is_attachable(
const struct gpart_fof_foreign *gp) {
return current_fof_attach_type & (1 << (gp->type + 1));
}
/**
* @brief Returns whether a foreign #gpart is of the 'linkable' kind.
*/
__attribute__((always_inline)) INLINE static int gpart_foreign_is_linkable(
const struct gpart_fof_foreign *gp) {
return current_fof_linking_type & (1 << (gp->type + 1));
}
/**
* @brief Returns whether a foreign #gpart is to be ignored by FOF.
*/
__attribute__((always_inline)) INLINE static int gpart_foreign_is_ignorable(
const struct gpart_fof_foreign *gp) {
return current_fof_ignore_type & (1 << (gp->type + 1));
}
/** /**
* @brief Finds the local root ID of the group a particle exists in. * @brief Finds the local root ID of the group a particle exists in.
* *
...@@ -1222,6 +1249,8 @@ static INLINE void add_foreign_link_to_list( ...@@ -1222,6 +1249,8 @@ static INLINE void add_foreign_link_to_list(
message("Re-allocating local group links from %d to %d elements.", message("Re-allocating local group links from %d to %d elements.",
*local_link_count, new_size); *local_link_count, new_size);
if (new_size < 0) error("Overflow in size of list of foreign links");
} }
/* Store the particle group properties for communication. */ /* Store the particle group properties for communication. */
...@@ -1250,7 +1279,7 @@ void fof_search_pair_cells_foreign( ...@@ -1250,7 +1279,7 @@ void fof_search_pair_cells_foreign(
const size_t count_i = ci->grav.count; const size_t count_i = ci->grav.count;
const size_t count_j = cj->grav.count; const size_t count_j = cj->grav.count;
const struct gpart *gparts_i = ci->grav.parts; const struct gpart *gparts_i = ci->grav.parts;
const struct gpart *gparts_j = cj->grav.parts; const struct gpart_fof_foreign *gparts_j = cj->grav.parts_fof_foreign;
/* Get local pointers */ /* Get local pointers */
const size_t *restrict group_index = props->group_index; const size_t *restrict group_index = props->group_index;
...@@ -1276,7 +1305,10 @@ void fof_search_pair_cells_foreign( ...@@ -1276,7 +1305,10 @@ void fof_search_pair_cells_foreign(
"cells."); "cells.");
if (!ci_local) { if (!ci_local) {
error("Cell ci, is not local."); error("Cell ci is not local!");
}
if (cj_local) {
error("Cell cj is local!");
} }
#endif #endif
...@@ -1322,16 +1354,16 @@ void fof_search_pair_cells_foreign( ...@@ -1322,16 +1354,16 @@ void fof_search_pair_cells_foreign(
for (size_t j = 0; j < count_j; j++) { for (size_t j = 0; j < count_j; j++) {
const struct gpart *pj = &gparts_j[j]; const struct gpart_fof_foreign *pj = &gparts_j[j];
/* Ignore inhibited particles */ /* Ignore inhibited particles */
if (pj->time_bin >= time_bin_inhibited) continue; if (pj->time_bin >= time_bin_inhibited) continue;
/* Check whether we ignore this particle type altogether */ /* Check whether we ignore this particle type altogether */
if (gpart_is_ignorable(pj)) continue; if (gpart_foreign_is_ignorable(pj)) continue;
/* Get the nature of the linking */ /* Get the nature of the linking */
const int is_link_j = gpart_is_linkable(pj); const int is_link_j = gpart_foreign_is_linkable(pj);
/* Only consider linkable<->linkable pairs */ /* Only consider linkable<->linkable pairs */
if (!(is_link_i && is_link_j)) continue; if (!(is_link_i && is_link_j)) continue;
...@@ -1560,20 +1592,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1560,20 +1592,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
const size_t count = c->grav.count; const size_t count = c->grav.count;
struct gpart *gparts = (struct gpart *)c->grav.parts; struct gpart *gparts = (struct gpart *)c->grav.parts;
/* Make a list of particle offsets into the global gparts array. */
size_t *const group_index = props->group_index;
#ifndef WITH_MPI
size_t *const index_offset = group_index + (ptrdiff_t)(gparts - space_gparts);
#endif
size_t *const attach_index = props->attach_index;
size_t *const attach_offset =
attach_index + (ptrdiff_t)(gparts - space_gparts);
char *const found_attach_index = props->found_attachable_link;
char *const found_attach_offset =
found_attach_index + (ptrdiff_t)(gparts - space_gparts);
/* Distances of particles in the global list */ /* Distances of particles in the global list */
float *const offset_dist = float *const offset_dist =
props->distance_to_link + (ptrdiff_t)(gparts - space_gparts); props->distance_to_link + (ptrdiff_t)(gparts - space_gparts);
...@@ -1603,14 +1621,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1603,14 +1621,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
const double piy = pi->x[1]; const double piy = pi->x[1];
const double piz = pi->x[2]; const double piz = pi->x[2];
/* Find the root of pi. */
#ifdef WITH_MPI
const size_t root_i = fof_find_global(
i + (ptrdiff_t)(gparts - space_gparts), group_index, nr_gparts);
#else
const size_t root_i = fof_find(index_offset[i], group_index);
#endif
/* Get the nature of the linking */ /* Get the nature of the linking */
const int is_link_i = gpart_is_linkable(pi); const int is_link_i = gpart_is_linkable(pi);
const int is_attach_i = gpart_is_attachable(pi); const int is_attach_i = gpart_is_attachable(pi);
...@@ -1648,14 +1658,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1648,14 +1658,6 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
error("Running FOF on an un-drifted particle!"); error("Running FOF on an un-drifted particle!");
#endif #endif
/* Find the root of pi. */
#ifdef WITH_MPI
const size_t root_j = fof_find_global(
j + (ptrdiff_t)(gparts - space_gparts), group_index, nr_gparts);
#else
const size_t root_j = fof_find(index_offset[j], group_index);
#endif
const double pjx = pj->x[0]; const double pjx = pj->x[0];
const double pjy = pj->x[1]; const double pjy = pj->x[1];
const double pjz = pj->x[2]; const double pjz = pj->x[2];
...@@ -1691,8 +1693,7 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1691,8 +1693,7 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
offset_dist[j] = dist; offset_dist[j] = dist;
/* Store the current best root */ /* Store the current best root */
attach_offset[j] = root_i; pj->fof_data.group_id = pi->fof_data.group_id;
found_attach_offset[j] = 1;
} }
} else if (is_link_j && is_attach_i) { } else if (is_link_j && is_attach_i) {
...@@ -1708,8 +1709,7 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1708,8 +1709,7 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
offset_dist[i] = dist; offset_dist[i] = dist;
/* Store the current best root */ /* Store the current best root */
attach_offset[i] = root_j; pi->fof_data.group_id = pj->fof_data.group_id;
found_attach_offset[i] = 1;
} }
} else { } else {
...@@ -1733,56 +1733,31 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2, ...@@ -1733,56 +1733,31 @@ void fof_attach_self_cell(const struct fof_props *props, const double l_x2,
* @param nr_gparts The number of #gpart in the local #space structure. * @param nr_gparts The number of #gpart in the local #space structure.
* @param ci The first #cell in which to perform FOF. * @param ci The first #cell in which to perform FOF.
* @param cj The second #cell in which to perform FOF. * @param cj The second #cell in which to perform FOF.
* @param ci_local Is the #cell ci on the local MPI rank?
* @param cj_local Is the #cell cj on the local MPI rank?
*/ */
void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], void fof_attach_pair_cells_both_local(const struct fof_props *props,
const double l_x2, const int periodic, const double dim[3], const double l_x2,
const struct gpart *const space_gparts, const int periodic,
const size_t nr_gparts, const struct gpart *const space_gparts,
const struct cell *restrict ci, const size_t nr_gparts,
const struct cell *restrict cj, const int ci_local, const struct cell *restrict ci,
const int cj_local) { const struct cell *restrict cj) {
#ifdef SWIFT_DEBUG_CHECKS
if (ci->nodeID != engine_rank) error("ci not local!");
if (cj->nodeID != engine_rank) error("cj not local!");
#endif
const size_t count_i = ci->grav.count; const size_t count_i = ci->grav.count;
const size_t count_j = cj->grav.count; const size_t count_j = cj->grav.count;
struct gpart *gparts_i = (struct gpart *)ci->grav.parts; struct gpart *gparts_i = (struct gpart *)ci->grav.parts;
struct gpart *gparts_j = (struct gpart *)cj->grav.parts; struct gpart *gparts_j = (struct gpart *)cj->grav.parts;
/* Index of particles in the global group list */
size_t *const group_index = props->group_index;
/* Make a list of particle offsets into the global gparts array. */
size_t *const index_offset_i =
group_index + (ptrdiff_t)(gparts_i - space_gparts);
size_t *const index_offset_j =
group_index + (ptrdiff_t)(gparts_j - space_gparts);
size_t *const attach_offset_i =
props->attach_index + (ptrdiff_t)(gparts_i - space_gparts);
size_t *const attach_offset_j =
props->attach_index + (ptrdiff_t)(gparts_j - space_gparts);
char *const found_attach_offset_i =
props->found_attachable_link + (ptrdiff_t)(gparts_i - space_gparts);
char *const found_attach_offset_j =
props->found_attachable_link + (ptrdiff_t)(gparts_j - space_gparts);
/* Distances of particles in the global list */ /* Distances of particles in the global list */
float *const offset_dist_i = float *const offset_dist_i =
props->distance_to_link + (ptrdiff_t)(gparts_i - space_gparts); props->distance_to_link + (ptrdiff_t)(gparts_i - space_gparts);
float *const offset_dist_j = float *const offset_dist_j =
props->distance_to_link + (ptrdiff_t)(gparts_j - space_gparts); props->distance_to_link + (ptrdiff_t)(gparts_j - space_gparts);
#ifdef SWIFT_DEBUG_CHECKS
if (index_offset_j > index_offset_i &&
(index_offset_j < index_offset_i + count_i))
error("Overlapping cells");
if (index_offset_i > index_offset_j &&
(index_offset_i < index_offset_j + count_j))
error("Overlapping cells");
#endif
/* Account for boundary conditions.*/ /* Account for boundary conditions.*/
double shift[3] = {0.0, 0.0, 0.0}; double shift[3] = {0.0, 0.0, 0.0};
...@@ -1819,19 +1794,6 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], ...@@ -1819,19 +1794,6 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
const double piy = pi->x[1] - shift[1]; const double piy = pi->x[1] - shift[1];
const double piz = pi->x[2] - shift[2]; const double piz = pi->x[2] - shift[2];
/* Find the root of pi. */
#ifdef WITH_MPI
size_t root_i;
if (ci_local) {
root_i = fof_find_global(index_offset_i[i] - node_offset, group_index,
nr_gparts);
} else {
root_i = pi->fof_data.group_id;
}
#else
const size_t root_i = fof_find(index_offset_i[i], group_index);
#endif
/* Get the nature of the linking */ /* Get the nature of the linking */
const int is_link_i = gpart_is_linkable(pi); const int is_link_i = gpart_is_linkable(pi);
const int is_attach_i = gpart_is_attachable(pi); const int is_attach_i = gpart_is_attachable(pi);
...@@ -1869,19 +1831,6 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], ...@@ -1869,19 +1831,6 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
error("Running FOF on an un-drifted particle!"); error("Running FOF on an un-drifted particle!");
#endif #endif
/* Find the root of pj. */
#ifdef WITH_MPI
size_t root_j;
if (cj_local) {
root_j = fof_find_global(index_offset_j[j] - node_offset, group_index,
nr_gparts);
} else {
root_j = pj->fof_data.group_id;
}
#else
const size_t root_j = fof_find(index_offset_j[j], group_index);
#endif
const double pjx = pj->x[0]; const double pjx = pj->x[0];
const double pjy = pj->x[1]; const double pjy = pj->x[1];
const double pjz = pj->x[2]; const double pjz = pj->x[2];
...@@ -1913,19 +1862,15 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], ...@@ -1913,19 +1862,15 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
* See whether it is closer and if so re-link. * See whether it is closer and if so re-link.
* This is safe to do as the attachables are never roots and * This is safe to do as the attachables are never roots and
* nothing is attached to them */ * nothing is attached to them */
if (cj_local) { const float dist = sqrtf(r2);
const float dist = sqrtf(r2);
if (dist < offset_dist_j[j]) { if (dist < offset_dist_j[j]) {
/* Store the new min dist */ /* Store the new min dist */
offset_dist_j[j] = dist; offset_dist_j[j] = dist;
/* Store the current best root */ /* Store the current best root */
attach_offset_j[j] = root_i; pj->fof_data.group_id = pi->fof_data.group_id;
found_attach_offset_j[j] = 1;
}
} }
} else if (is_link_j && is_attach_i) { } else if (is_link_j && is_attach_i) {
...@@ -1934,19 +1879,15 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], ...@@ -1934,19 +1879,15 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
* See whether it is closer and if so re-link. * See whether it is closer and if so re-link.
* This is safe to do as the attachables are never roots and * This is safe to do as the attachables are never roots and
* nothing is attached to them */ * nothing is attached to them */
if (ci_local) { const float dist = sqrtf(r2);
const float dist = sqrtf(r2);
if (dist < offset_dist_i[i]) { if (dist < offset_dist_i[i]) {
/* Store the new min dist */ /* Store the new min dist */
offset_dist_i[i] = dist; offset_dist_i[i] = dist;
/* Store the current best root */ /* Store the current best root */
attach_offset_i[i] = root_j; pi->fof_data.group_id = pj->fof_data.group_id;
found_attach_offset_i[i] = 1;
}
} }
} else { } else {
...@@ -1959,762 +1900,486 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3], ...@@ -1959,762 +1900,486 @@ void fof_attach_pair_cells(const struct fof_props *props, const double dim[3],
} }
} }
/** void fof_attach_pair_cells_ci_local(const struct fof_props *props,
* @brief Recursively perform a union-find attaching between two cells. const double dim[3], const double l_x2,
* const int periodic,
* If cells are more distant than the linking length, we abort early. const struct gpart *const space_gparts,
* const size_t nr_gparts,
* @param props The properties fof the FOF scheme. const struct cell *restrict ci,
* @param dim The dimension of the space. const struct cell *restrict cj) {
* @param attach_r2 the square of the FOF linking length.
* @param periodic Are we using periodic BCs?
* @param space_gparts The start of the #gpart array in the #space structure.
* @param nr_gparts The number of #gpart in the local #space structure.
* @param ci The first #cell in which to perform FOF.
* @param cj The second #cell in which to perform FOF.
* @param ci_local Is the #cell ci on the local MPI rank?
* @param cj_local Is the #cell cj on the local MPI rank?
*/
void rec_fof_attach_pair(const struct fof_props *props, const double dim[3],
const double attach_r2, const int periodic,
const struct gpart *const space_gparts,
const size_t nr_gparts, struct cell *restrict ci,
struct cell *restrict cj, const int ci_local,
const int cj_local) {
/* Find the shortest distance between cells, remembering to account for
* boundary conditions. */
const double r2 = cell_min_dist(ci, cj, dim);
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
if (ci == cj) error("Pair FOF called on same cell!!!"); if (ci->nodeID != engine_rank) error("ci not local!");
if (cj->nodeID == engine_rank) error("cj local!");
#endif #endif
/* Return if cells are out of range of each other. */ const size_t count_i = ci->grav.count;
if (r2 > attach_r2) return; const size_t count_j = cj->grav.count;
struct gpart *gparts_i = ci->grav.parts;
struct gpart_fof_foreign *gparts_j = cj->grav.parts_fof_foreign;
/* Recurse on both cells if they are both split. */ /* Distances of particles in the global list */
if (ci->split && cj->split) { float *const offset_dist_i =
for (int k = 0; k < 8; k++) { props->distance_to_link + (ptrdiff_t)(gparts_i - space_gparts);
if (ci->progeny[k] != NULL) {
for (int l = 0; l < 8; l++) /* Account for boundary conditions.*/
if (cj->progeny[l] != NULL) double shift[3] = {0.0, 0.0, 0.0};
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, ci->progeny[k], cj->progeny[l], /* Get the relative distance between the pairs, wrapping. */
ci_local, cj_local); double diff[3];
} for (int k = 0; k < 3; k++) {
} diff[k] = cj->loc[k] - ci->loc[k];
} else if (ci->split) { if (periodic && diff[k] < -dim[k] * 0.5)
for (int k = 0; k < 8; k++) { shift[k] = dim[k];
if (ci->progeny[k] != NULL) else if (periodic && diff[k] > dim[k] * 0.5)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts, shift[k] = -dim[k];
nr_gparts, ci->progeny[k], cj, ci_local, cj_local); else
} shift[k] = 0.0;
} else if (cj->split) { diff[k] += shift[k];
for (int k = 0; k < 8; k++) {
if (cj->progeny[k] != NULL)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, ci, cj->progeny[k], ci_local, cj_local);
}
} else {
/* Perform FOF attach between pairs of cells that are within the linking
* length and not the same cell. */
fof_attach_pair_cells(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, ci, cj, ci_local, cj_local);
} }
}
/** /* Loop over particles and find which particles belong in the same group. */
* @brief Recursively perform a the attaching operation on a cell. for (size_t i = 0; i < count_i; i++) {
*
* @param props The properties fof the FOF scheme.
* @param dim The dimension of the space.
* @param attach_r2 the square of the FOF linking length.
* @param periodic Are we using periodic BCs?
* @param space_gparts The start of the #gpart array in the #space structure.
* @param nr_gparts The number of #gpart in the local #space structure.
* @param c The #cell in which to perform FOF.
*/
void rec_fof_attach_self(const struct fof_props *props, const double dim[3],
const double attach_r2, const int periodic,
const struct gpart *const space_gparts,
const size_t nr_gparts, struct cell *c) {
/* Recurse? */ struct gpart *restrict pi = &gparts_i[i];
if (c->split) {
/* Loop over all progeny. Perform pair and self recursion on progenies.*/ /* Ignore inhibited particles */
for (int k = 0; k < 8; k++) { if (pi->time_bin >= time_bin_inhibited) continue;
if (c->progeny[k] != NULL) {
rec_fof_attach_self(props, dim, attach_r2, periodic, space_gparts, /* Check whether we ignore this particle type altogether */
nr_gparts, c->progeny[k]); if (gpart_is_ignorable(pi)) continue;
for (int l = k + 1; l < 8; l++) #ifdef SWIFT_DEBUG_CHECKS
if (c->progeny[l] != NULL) if (pi->ti_drift != ti_current)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts, error("Running FOF on an un-drifted particle!");
nr_gparts, c->progeny[k], c->progeny[l], #endif
/*ci_local=*/1,
/*cj_local=*/1);
}
}
} else {
/* Otherwise, compute self-interaction. */ const double pix = pi->x[0] - shift[0];
fof_attach_self_cell(props, attach_r2, space_gparts, nr_gparts, c); const double piy = pi->x[1] - shift[1];
} const double piz = pi->x[2] - shift[2];
}
/* Mapper function to atomically update the group size array. */ /* Get the nature of the linking */
void fof_update_group_size_mapper(hashmap_key_t key, hashmap_value_t *value, const int is_link_i = gpart_is_linkable(pi);
void *data) { const int is_attach_i = gpart_is_attachable(pi);
size_t *group_size = (size_t *)data; #ifdef SWIFT_DEBUG_CHECKS
if (is_link_i && is_attach_i)
error("Particle cannot be both linkable and attachable!");
#endif
/* Use key to index into group size array. */ for (size_t j = 0; j < count_j; j++) {
atomic_add(&group_size[key], value->value_st);
}
/** struct gpart_fof_foreign *restrict pj = &gparts_j[j];
* @brief Mapper function to calculate the group sizes.
*
* @param map_data An array of #gpart%s.
* @param num_elements Chunk size.
* @param extra_data Pointer to a #space.
*/
void fof_calc_group_size_mapper(void *map_data, int num_elements,
void *extra_data) {
/* Retrieve mapped data. */ /* Ignore inhibited particles */
struct space *s = (struct space *)extra_data; if (pj->time_bin >= time_bin_inhibited) continue;
struct gpart *gparts = (struct gpart *)map_data;
size_t *restrict group_index = s->e->fof_properties->group_index;
size_t *restrict group_size = s->e->fof_properties->group_size;
/* Offset into gparts array. */ /* Check whether we ignore this particle type altogether */
const ptrdiff_t gparts_offset = (ptrdiff_t)(gparts - s->gparts); if (gpart_foreign_is_ignorable(pj)) continue;
size_t *const group_index_offset = group_index + gparts_offset;
/* Create hash table. */ /* Get the nature of the linking */
hashmap_t map; const int is_link_j = gpart_foreign_is_linkable(pj);
hashmap_init(&map); const int is_attach_j = gpart_foreign_is_attachable(pj);
for (int ind = 0; ind < num_elements; ind++) { #ifdef SWIFT_DEBUG_CHECKS
if (is_link_j && is_attach_j)
error("Particle cannot be both linkable and attachable!");
#endif
const hashmap_key_t root = /* We only want link<->attach pairs */
(hashmap_key_t)fof_find(group_index_offset[ind], group_index); if (is_attach_i && is_attach_j) continue;
const size_t gpart_index = gparts_offset + ind; if (is_link_i && is_link_j) continue;
/* Only add particles which aren't the root of a group. Stops groups of size #ifdef SWIFT_DEBUG_CHECKS
* 1 being added to the hash table. */ if (pj->ti_drift != ti_current)
if (root != gpart_index) { error("Running FOF on an un-drifted particle!");
hashmap_value_t *size = hashmap_get(&map, root); #endif
if (size != NULL)
(*size).value_st++;
else
error("Couldn't find key (%zu) or create new one.", root);
}
}
/* Update the group size array. */ const double pjx = pj->x[0];
if (map.size > 0) const double pjy = pj->x[1];
hashmap_iterate(&map, fof_update_group_size_mapper, group_size); const double pjz = pj->x[2];
hashmap_free(&map); /* Compute pairwise distance (periodic BCs were accounted
} for by the shift vector) */
float dx[3], r2 = 0.0f;
dx[0] = pix - pjx;
dx[1] = piy - pjy;
dx[2] = piz - pjz;
/* Mapper function to atomically update the group mass array. */ for (int k = 0; k < 3; k++) r2 += dx[k] * dx[k];
static INLINE void fof_update_group_mass_iterator(hashmap_key_t key,
hashmap_value_t *value,
void *data) {
double *group_mass = (double *)data; /* Hit or miss? */
if (r2 < l_x2) {
/* Use key to index into group mass array. */ /* Now that we are within the linking length,
atomic_add_d(&group_mass[key], value->value_dbl); * decide what to do based on linking types */
}
/* Mapper function to atomically update the group size array. */ if (is_link_i && is_link_j) {
static INLINE void fof_update_group_size_iterator(hashmap_key_t key,
hashmap_value_t *value,
void *data) {
long long *group_size = (long long *)data;
/* Use key to index into group mass array. */ #ifdef SWIFT_DEBUG_CHECKS
atomic_add(&group_size[key], value->value_st); error("Fundamental logic error!");
} #endif
/** } else if (is_link_i && is_attach_j) {
* @brief Mapper function to calculate the group masses.
*
* @param map_data An array of #gpart%s.
* @param num_elements Chunk size.
* @param extra_data Pointer to a #space.
*/
void fof_calc_group_mass_mapper(void *map_data, int num_elements,
void *extra_data) {
/* Retrieve mapped data. */ /* Nothing to do here. The reverse action will be done
struct space *s = (struct space *)extra_data; in the converse call on the other node. */
struct gpart *gparts = (struct gpart *)map_data;
double *group_mass = s->e->fof_properties->group_mass;
long long *group_size = s->e->fof_properties->final_group_size;
const size_t group_id_default = s->e->fof_properties->group_id_default;
const size_t group_id_offset = s->e->fof_properties->group_id_offset;
/* Create hash table. */ } else if (is_link_j && is_attach_i) {
hashmap_t map;
hashmap_init(&map);
/* Loop over particles and increment the group mass for groups above /* We got a linkable and an attachable.
* min_group_size. */ * See whether it is closer and if so re-link.
for (int ind = 0; ind < num_elements; ind++) { * This is safe to do as the attachables are never roots and
* nothing is attached to them */
const float dist = sqrtf(r2);
/* Only check groups above the minimum size. */ if (dist < offset_dist_i[i]) {
if (gparts[ind].fof_data.group_id != group_id_default) {
hashmap_key_t index = gparts[ind].fof_data.group_id - group_id_offset; /* Store the new min dist */
hashmap_value_t *data = hashmap_get(&map, index); offset_dist_i[i] = dist;
/* Update group mass */ /* Store the current best root */
if (data != NULL) { pi->fof_data.group_id = pj->fof_data.group_id;
(*data).value_dbl += gparts[ind].mass; }
(*data).value_st += 1; } else {
} else #ifdef SWIFT_DEBUG_CHECKS
error("Couldn't find key (%zu) or create new one.", index); error("Fundamental logic error!");
#endif
}
}
} }
} }
/* Update the group mass array. */
if (map.size > 0)
hashmap_iterate(&map, fof_update_group_mass_iterator, group_mass);
if (map.size > 0)
hashmap_iterate(&map, fof_update_group_size_iterator, group_size);
hashmap_free(&map);
} }
#ifdef WITH_MPI void fof_attach_pair_cells_cj_local(const struct fof_props *props,
/* Mapper function to unpack hash table into array. */ const double dim[3], const double l_x2,
void fof_unpack_group_mass_mapper(hashmap_key_t key, hashmap_value_t *value, const int periodic,
void *data) { const struct gpart *const space_gparts,
const size_t nr_gparts,
struct fof_mass_send_hashmap *fof_mass_send = const struct cell *restrict ci,
(struct fof_mass_send_hashmap *)data; const struct cell *restrict cj) {
struct fof_final_mass *mass_send = fof_mass_send->mass_send;
size_t *nsend = &fof_mass_send->nsend;
/* Store elements from hash table in array. */
mass_send[*nsend].global_root = key;
mass_send[*nsend].group_mass = value->value_dbl;
mass_send[*nsend].final_group_size = value->value_ll;
mass_send[*nsend].first_position[0] = value->value_array2_dbl[0];
mass_send[*nsend].first_position[1] = value->value_array2_dbl[1];
mass_send[*nsend].first_position[2] = value->value_array2_dbl[2];
mass_send[*nsend].centre_of_mass[0] = value->value_array_dbl[0];
mass_send[*nsend].centre_of_mass[1] = value->value_array_dbl[1];
mass_send[*nsend].centre_of_mass[2] = value->value_array_dbl[2];
mass_send[*nsend].max_part_density_index = value->value_st;
mass_send[*nsend].max_part_density = value->value_flt;
(*nsend)++;
}
#endif /* WITH_MPI */ #ifdef SWIFT_DEBUG_CHECKS
if (ci->nodeID == engine_rank) error("ci local!");
if (cj->nodeID != engine_rank) error("cj not local!");
#endif
/** const size_t count_i = ci->grav.count;
* @brief Calculates the total mass and CoM of each group above min_group_size const size_t count_j = cj->grav.count;
* and finds the densest particle for black hole seeding. struct gpart_fof_foreign *gparts_i = ci->grav.parts_fof_foreign;
*/ struct gpart *gparts_j = cj->grav.parts;
void fof_calc_group_mass(struct fof_props *props, const struct space *s,
const int seed_black_holes,
const size_t num_groups_local,
const size_t num_groups_prev,
size_t *restrict num_on_node,
size_t *restrict first_on_node,
double *restrict group_mass) {
const size_t nr_gparts = s->nr_gparts; /* Distances of particles in the global list */
struct gpart *gparts = s->gparts; float *const offset_dist_j =
const struct part *parts = s->parts; props->distance_to_link + (ptrdiff_t)(gparts_j - space_gparts);
const size_t group_id_offset = props->group_id_offset;
const size_t group_id_default = props->group_id_default;
const double seed_halo_mass = props->seed_halo_mass;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
#ifdef WITH_MPI /* Account for boundary conditions.*/
size_t *group_index = props->group_index; double shift[3] = {0.0, 0.0, 0.0};
const int nr_nodes = s->e->nr_nodes;
/* Direct pointers to the arrays */ /* Get the relative distance between the pairs, wrapping. */
long long *max_part_density_index = props->max_part_density_index; double diff[3];
float *max_part_density = props->max_part_density; for (int k = 0; k < 3; k++) {
double *centre_of_mass = props->group_centre_of_mass; diff[k] = cj->loc[k] - ci->loc[k];
double *first_position = props->group_first_position; if (periodic && diff[k] < -dim[k] * 0.5)
long long *final_group_size = props->final_group_size; shift[k] = dim[k];
else if (periodic && diff[k] > dim[k] * 0.5)
shift[k] = -dim[k];
else
shift[k] = 0.0;
diff[k] += shift[k];
}
/* Start the hash map */ /* Loop over particles and find which particles belong in the same group. */
hashmap_t map; for (size_t i = 0; i < count_i; i++) {
hashmap_init(&map);
/* Collect information about the local particles and update the local AND struct gpart_fof_foreign *restrict pi = &gparts_i[i];
* foreign group fragments */
for (size_t i = 0; i < nr_gparts; i++) {
/* Ignore inhibited particles */ /* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue; if (pi->time_bin >= time_bin_inhibited) continue;
/* Check whether we ignore this particle type altogether */ /* Check whether we ignore this particle type altogether */
if (gpart_is_ignorable(&gparts[i])) continue; if (gpart_foreign_is_ignorable(pi)) continue;
/* Check if the particle is in a group above the threshold. */
if (gparts[i].fof_data.group_id != group_id_default) {
const size_t root = fof_find_global(i, group_index, nr_gparts);
if (is_local(root, nr_gparts)) {
/* The root is local */ #ifdef SWIFT_DEBUG_CHECKS
if (pi->ti_drift != ti_current)
const size_t index = error("Running FOF on an un-drifted particle!");
gparts[i].fof_data.group_id - group_id_offset - num_groups_prev; #endif
/* Update group mass */
group_mass[index] += gparts[i].mass;
/* Update group size */
final_group_size[index]++;
} else {
/* The root is *not* local */
/* Get the root in the foreign hashmap (create if necessary) */
hashmap_value_t *const data = hashmap_get(&map, (hashmap_key_t)root);
if (data == NULL)
error("Couldn't find key (%zu) or create new one.", root);
/* Compute the centre of mass */
const double mass = gparts[i].mass;
double x[3] = {gparts[i].x[0], gparts[i].x[1], gparts[i].x[2]};
/* Add mass fragments of groups */
data->value_dbl += mass;
/* Increase fragment size */
data->value_ll++;
/* Record the first particle of this fragment that we encounter so we
* we can use it as reference frame for the centre of mass calculation
*/
if (data->value_array2_dbl[0] == (double)(-FLT_MAX)) {
data->value_array2_dbl[0] = gparts[i].x[0];
data->value_array2_dbl[1] = gparts[i].x[1];
data->value_array2_dbl[2] = gparts[i].x[2];
}
if (periodic) {
x[0] = nearest(x[0] - data->value_array2_dbl[0], dim[0]);
x[1] = nearest(x[1] - data->value_array2_dbl[1], dim[1]);
x[2] = nearest(x[2] - data->value_array2_dbl[2], dim[2]);
}
data->value_array_dbl[0] += mass * x[0];
data->value_array_dbl[1] += mass * x[1];
data->value_array_dbl[2] += mass * x[2];
/* Also accumulate the densest gas particle and its index */
if (gparts[i].type == swift_type_gas &&
data->value_st != fof_halo_has_black_hole) {
const size_t gas_index = -gparts[i].id_or_neg_offset;
const float rho_com = hydro_get_comoving_density(&parts[gas_index]);
/* Update index if a denser gas particle is found. */
if (rho_com > data->value_flt) {
data->value_flt = rho_com;
data->value_st = gas_index;
}
} else if (gparts[i].type == swift_type_black_hole) {
/* If there is already a black hole in the fragment we don't need to
create a new one. */
data->value_st = fof_halo_has_black_hole;
data->value_flt = 0.f;
}
} /* Foreign root */
} /* Particle is in a group */
} /* Loop over particles */
size_t nsend = map.size;
struct fof_mass_send_hashmap hashmap_mass_send = {NULL, 0};
/* Allocate and initialise a mass array. */
if (posix_memalign((void **)&hashmap_mass_send.mass_send, 32,
nsend * sizeof(struct fof_final_mass)) != 0)
error("Failed to allocate list of group masses for FOF search.");
hashmap_mass_send.nsend = 0;
struct fof_final_mass *fof_mass_send = hashmap_mass_send.mass_send;
/* Unpack mass fragments and roots from hash table. */ const double pix = pi->x[0] - shift[0];
if (map.size > 0) const double piy = pi->x[1] - shift[1];
hashmap_iterate(&map, fof_unpack_group_mass_mapper, &hashmap_mass_send); const double piz = pi->x[2] - shift[2];
nsend = hashmap_mass_send.nsend; /* Get the nature of the linking */
const int is_link_i = gpart_foreign_is_linkable(pi);
const int is_attach_i = gpart_foreign_is_attachable(pi);
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
if (nsend != map.size) if (is_link_i && is_attach_i)
error("No. of mass fragments to send != elements in hash table."); error("Particle cannot be both linkable and attachable!");
#endif #endif
hashmap_free(&map); for (size_t j = 0; j < count_j; j++) {
/* Sort by global root - this puts the groups in order of which node they're struct gpart *restrict pj = &gparts_j[j];
* stored on */
qsort(fof_mass_send, nsend, sizeof(struct fof_final_mass),
compare_fof_final_mass_global_root);
/* Determine how many entries go to each node */ /* Ignore inhibited particles */
int *sendcount = (int *)calloc(nr_nodes, sizeof(int)); if (pj->time_bin >= time_bin_inhibited) continue;
int dest = 0;
for (size_t i = 0; i < nsend; i += 1) {
while ((fof_mass_send[i].global_root >=
first_on_node[dest] + num_on_node[dest]) ||
(num_on_node[dest] == 0))
dest += 1;
if (dest >= nr_nodes) error("Node index out of range!");
sendcount[dest] += 1;
}
int *recvcount = NULL, *sendoffset = NULL, *recvoffset = NULL; /* Check whether we ignore this particle type altogether */
size_t nrecv = 0; if (gpart_is_ignorable(pj)) continue;
fof_compute_send_recv_offsets(nr_nodes, sendcount, &recvcount, &sendoffset, /* Get the nature of the linking */
&recvoffset, &nrecv); const int is_link_j = gpart_is_linkable(pj);
const int is_attach_j = gpart_is_attachable(pj);
struct fof_final_mass *fof_mass_recv = #ifdef SWIFT_DEBUG_CHECKS
(struct fof_final_mass *)malloc(nrecv * sizeof(struct fof_final_mass)); if (is_link_j && is_attach_j)
error("Particle cannot be both linkable and attachable!");
#endif
/* Exchange group mass */ /* We only want link<->attach pairs */
MPI_Alltoallv(fof_mass_send, sendcount, sendoffset, fof_final_mass_type, if (is_attach_i && is_attach_j) continue;
fof_mass_recv, recvcount, recvoffset, fof_final_mass_type, if (is_link_i && is_link_j) continue;
MPI_COMM_WORLD);
/* For each received global root, look up the group ID we assigned and
* increment the group mass */
for (size_t i = 0; i < nrecv; i++) {
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
if ((fof_mass_recv[i].global_root < node_offset) || if (pj->ti_drift != ti_current)
(fof_mass_recv[i].global_root >= node_offset + nr_gparts)) { error("Running FOF on an un-drifted particle!");
error("Received global root index out of range!");
}
#endif #endif
const size_t local_root_index = fof_mass_recv[i].global_root - node_offset;
const size_t local_group_offset = group_id_offset + num_groups_prev;
const size_t index =
gparts[local_root_index].fof_data.group_id - local_group_offset;
group_mass[index] += fof_mass_recv[i].group_mass;
final_group_size[index] += fof_mass_recv[i].final_group_size;
}
/* Loop over particles, densest particle in each *local* group. const double pjx = pj->x[0];
* We can do this now as we eventually have the total group mass */ const double pjy = pj->x[1];
for (size_t i = 0; i < nr_gparts; i++) { const double pjz = pj->x[2];
/* Ignore inhibited particles */ /* Compute pairwise distance (periodic BCs were accounted
if (gparts[i].time_bin >= time_bin_inhibited) continue; for by the shift vector) */
float dx[3], r2 = 0.0f;
dx[0] = pix - pjx;
dx[1] = piy - pjy;
dx[2] = piz - pjz;
/* Check whether we ignore this particle type altogether */ for (int k = 0; k < 3; k++) r2 += dx[k] * dx[k];
if (current_fof_ignore_type & (1 << (gparts[i].type + 1))) continue;
/* Only check groups above the minimum mass threshold. */ /* Hit or miss? */
if (gparts[i].fof_data.group_id != group_id_default) { if (r2 < l_x2) {
const size_t root = fof_find_global(i, group_index, nr_gparts); /* Now that we are within the linking length,
* decide what to do based on linking types */
if (is_local(root, nr_gparts)) { if (is_link_i && is_link_j) {
const size_t index = #ifdef SWIFT_DEBUG_CHECKS
gparts[i].fof_data.group_id - group_id_offset - num_groups_prev; error("Fundamental logic error!");
#endif
/* Compute the centre of mass */ } else if (is_link_i && is_attach_j) {
const double mass = gparts[i].mass;
double x[3] = {gparts[i].x[0], gparts[i].x[1], gparts[i].x[2]};
/* Record the first particle of this group that we encounter so we /* We got a linkable and an attachable.
* can use it as reference frame for the centre of mass calculation */ * See whether it is closer and if so re-link.
if (first_position[index * 3 + 0] == (double)(-FLT_MAX)) { * This is safe to do as the attachables are never roots and
first_position[index * 3 + 0] = x[0]; * nothing is attached to them */
first_position[index * 3 + 1] = x[1];
first_position[index * 3 + 2] = x[2];
}
if (periodic) { const float dist = sqrtf(r2);
x[0] = nearest(x[0] - first_position[index * 3 + 0], dim[0]);
x[1] = nearest(x[1] - first_position[index * 3 + 1], dim[1]);
x[2] = nearest(x[2] - first_position[index * 3 + 2], dim[2]);
}
centre_of_mass[index * 3 + 0] += mass * x[0]; if (dist < offset_dist_j[j]) {
centre_of_mass[index * 3 + 1] += mass * x[1];
centre_of_mass[index * 3 + 2] += mass * x[2];
/* Check haloes above the seeding threshold */ /* Store the new min dist */
if (group_mass[index] > seed_halo_mass) { offset_dist_j[j] = dist;
/* Find the densest gas particle. /* Store the current best root */
* Account for groups that already have a black hole and groups that pj->fof_data.group_id = pi->fof_data.group_id;
* contain no gas. */ }
if (gparts[i].type == swift_type_gas &&
max_part_density_index[index] != fof_halo_has_black_hole) {
const size_t gas_index = -gparts[i].id_or_neg_offset; } else if (is_link_j && is_attach_i) {
const float rho_com = hydro_get_comoving_density(&parts[gas_index]);
/* Update index if a denser gas particle is found. */ /* Nothing to do here. The reverse action will be done
if (rho_com > max_part_density[index]) { in the converse call on the other node. */
max_part_density_index[index] = gas_index;
max_part_density[index] = rho_com;
}
}
/* If there is already a black hole in the group we don't need to
create a new one. */
else if (gparts[i].type == swift_type_black_hole) {
max_part_density_index[index] = fof_halo_has_black_hole;
}
} else { } else {
max_part_density_index[index] = fof_halo_has_too_low_mass; #ifdef SWIFT_DEBUG_CHECKS
error("Fundamental logic error!");
#endif
} }
} }
} }
} }
}
/* For each received global root, look up the group ID we assigned and find /**
* the global maximum gas density */ * @brief Recursively perform a union-find attaching between two cells.
for (size_t i = 0; i < nrecv; i++) { *
* If cells are more distant than the linking length, we abort early.
const size_t local_root_index = fof_mass_recv[i].global_root - node_offset; *
const size_t local_group_offset = group_id_offset + num_groups_prev; * @param props The properties fof the FOF scheme.
const size_t index = * @param dim The dimension of the space.
gparts[local_root_index].fof_data.group_id - local_group_offset; * @param attach_r2 the square of the FOF linking length.
* @param periodic Are we using periodic BCs?
double fragment_mass = fof_mass_recv[i].group_mass; * @param space_gparts The start of the #gpart array in the #space structure.
double fragment_centre_of_mass[3] = { * @param nr_gparts The number of #gpart in the local #space structure.
fof_mass_recv[i].centre_of_mass[0] / fof_mass_recv[i].group_mass, * @param ci The first #cell in which to perform FOF.
fof_mass_recv[i].centre_of_mass[1] / fof_mass_recv[i].group_mass, * @param cj The second #cell in which to perform FOF.
fof_mass_recv[i].centre_of_mass[2] / fof_mass_recv[i].group_mass}; * @param ci_local Is the #cell ci on the local MPI rank?
fragment_centre_of_mass[0] += fof_mass_recv[i].first_position[0]; * @param cj_local Is the #cell cj on the local MPI rank?
fragment_centre_of_mass[1] += fof_mass_recv[i].first_position[1]; */
fragment_centre_of_mass[2] += fof_mass_recv[i].first_position[2]; void rec_fof_attach_pair(const struct fof_props *props, const double dim[3],
const double attach_r2, const int periodic,
const struct gpart *const space_gparts,
const size_t nr_gparts, struct cell *restrict ci,
struct cell *restrict cj, const int ci_local,
const int cj_local) {
if (periodic) { /* Find the shortest distance between cells, remembering to account for
fragment_centre_of_mass[0] = nearest( * boundary conditions. */
fragment_centre_of_mass[0] - first_position[3 * index + 0], dim[0]); const double r2 = cell_min_dist(ci, cj, dim);
fragment_centre_of_mass[1] = nearest(
fragment_centre_of_mass[1] - first_position[3 * index + 1], dim[1]);
fragment_centre_of_mass[2] = nearest(
fragment_centre_of_mass[2] - first_position[3 * index + 2], dim[2]);
}
centre_of_mass[index * 3 + 0] += fragment_mass * fragment_centre_of_mass[0]; #ifdef SWIFT_DEBUG_CHECKS
centre_of_mass[index * 3 + 1] += fragment_mass * fragment_centre_of_mass[1]; if (ci == cj) error("Pair FOF called on same cell!!!");
centre_of_mass[index * 3 + 2] += fragment_mass * fragment_centre_of_mass[2]; #endif
/* Only seed groups above the mass threshold. */ /* Return if cells are out of range of each other. */
if (group_mass[index] > seed_halo_mass) { if (r2 > attach_r2) return;
/* Only check groups that don't already contain a black hole. */ /* Recurse on both cells if they are both split. */
if (max_part_density_index[index] != fof_halo_has_black_hole) { if (ci->split && cj->split) {
for (int k = 0; k < 8; k++) {
if (ci->progeny[k] != NULL) {
/* Find the densest particle in each group using the densest particle for (int l = 0; l < 8; l++)
* from each group fragment. */ if (cj->progeny[l] != NULL)
if (fof_mass_recv[i].max_part_density > max_part_density[index]) { rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
max_part_density[index] = fof_mass_recv[i].max_part_density; nr_gparts, ci->progeny[k], cj->progeny[l],
max_part_density_index[index] = ci_local, cj_local);
fof_mass_recv[i].max_part_density_index;
}
}
/* If there is already a black hole in the group we don't need to create a
new one. */
else if (fof_mass_recv[i].max_part_density_index ==
fof_halo_has_black_hole) {
max_part_density_index[index] = fof_halo_has_black_hole;
} }
} else {
max_part_density_index[index] = fof_halo_has_too_low_mass;
} }
} } else if (ci->split) {
for (int k = 0; k < 8; k++) {
if (ci->progeny[k] != NULL)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, ci->progeny[k], cj, ci_local, cj_local);
}
} else if (cj->split) {
for (int k = 0; k < 8; k++) {
if (cj->progeny[k] != NULL)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, ci, cj->progeny[k], ci_local, cj_local);
}
} else {
/* Perform FOF attach between pairs of cells that are within the linking
* length and not the same cell. */
if (ci_local && cj_local) {
/* For each received global root, look up the group ID we assigned and send fof_attach_pair_cells_both_local(props, dim, attach_r2, periodic,
* the global maximum gas density index back */ space_gparts, nr_gparts, ci, cj);
for (size_t i = 0; i < nrecv; i++) { } else if (ci_local && !cj_local) {
fof_attach_pair_cells_ci_local(props, dim, attach_r2, periodic,
space_gparts, nr_gparts, ci, cj);
#ifdef SWIFT_DEBUG_CHECKS } else if (!ci_local && cj_local) {
if ((fof_mass_recv[i].global_root < node_offset) || fof_attach_pair_cells_cj_local(props, dim, attach_r2, periodic,
(fof_mass_recv[i].global_root >= node_offset + nr_gparts)) { space_gparts, nr_gparts, ci, cj);
error("Received global root index out of range!");
}
#endif
const size_t local_root_index = fof_mass_recv[i].global_root - node_offset;
const size_t local_group_offset = group_id_offset + num_groups_prev;
const size_t index =
gparts[local_root_index].fof_data.group_id - local_group_offset;
/* If the densest particle found locally is not the global max, make sure we
* don't seed two black holes. */
if (max_part_density_index[index] ==
fof_mass_recv[i].max_part_density_index) {
/* If the local index has been set to a foreign index then we don't need
* to seed a black hole locally. */
max_part_density_index[index] = fof_halo_has_black_hole;
} else { } else {
/* The densest particle is on the same node as the global root so we don't error("Error in the recursion logic");
need to seed a black hole on the other node. */
fof_mass_recv[i].max_part_density_index = fof_halo_has_black_hole;
} }
} }
}
/* Send the result back */ /**
MPI_Alltoallv(fof_mass_recv, recvcount, recvoffset, fof_final_mass_type, * @brief Recursively perform a the attaching operation on a cell.
fof_mass_send, sendcount, sendoffset, fof_final_mass_type, *
MPI_COMM_WORLD); * @param props The properties fof the FOF scheme.
* @param dim The dimension of the space.
int extra_seed_count = 0; * @param attach_r2 the square of the FOF linking length.
size_t density_index_size = num_groups_local; * @param periodic Are we using periodic BCs?
* @param space_gparts The start of the #gpart array in the #space structure.
/* Add the index of the densest particle to the local list if the global root * @param nr_gparts The number of #gpart in the local #space structure.
* is not on this node. */ * @param c The #cell in which to perform FOF.
for (size_t i = 0; i < nsend; i++) { */
void rec_fof_attach_self(const struct fof_props *props, const double dim[3],
/* Only add the index if: const double attach_r2, const int periodic,
* 1) there is not already a black hole in the group const struct gpart *const space_gparts,
* AND const size_t nr_gparts, struct cell *c) {
* 2) there is gas in the group. */
if (fof_mass_send[i].max_part_density_index >= 0) {
/* Re-allocate the list if it's needed. */
if (num_groups_local + extra_seed_count >= density_index_size) {
const size_t new_size = 2 * density_index_size;
max_part_density_index = (long long *)realloc( /* Recurse? */
max_part_density_index, new_size * sizeof(long long)); if (c->split) {
message( /* Loop over all progeny. Perform pair and self recursion on progenies.*/
"Re-allocating max_part_density_index from %zu to %zu elements.", for (int k = 0; k < 8; k++) {
density_index_size, new_size); if (c->progeny[k] != NULL) {
density_index_size = new_size; rec_fof_attach_self(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, c->progeny[k]);
props->max_part_density_index = max_part_density_index; for (int l = k + 1; l < 8; l++)
if (c->progeny[l] != NULL)
rec_fof_attach_pair(props, dim, attach_r2, periodic, space_gparts,
nr_gparts, c->progeny[k], c->progeny[l],
/*ci_local=*/1,
/*cj_local=*/1);
} }
/* Add particle index onto the end of the array. */
max_part_density_index[num_groups_local + extra_seed_count] =
fof_mass_send[i].max_part_density_index;
extra_seed_count++;
} }
} } else {
props->extra_bh_seed_count = extra_seed_count;
free(sendcount);
free(recvcount);
free(sendoffset);
free(recvoffset);
free(fof_mass_send);
free(fof_mass_recv);
#else
/* Increment the group mass for groups above min_group_size. */
threadpool_map(&s->e->threadpool, fof_calc_group_mass_mapper, gparts,
nr_gparts, sizeof(struct gpart), threadpool_auto_chunk_size,
(struct space *)s);
/* Direct pointers to the arrays */
long long *max_part_density_index = props->max_part_density_index;
float *max_part_density = props->max_part_density;
double *centre_of_mass = props->group_centre_of_mass;
double *first_position = props->group_first_position;
/* Loop over particles, compute CoM and find the densest particle in each
* group. */
/* JSW TODO: Parallelise with threadpool*/
for (size_t i = 0; i < nr_gparts; i++) {
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
/* Check whether we ignore this particle type altogether */ /* Otherwise, compute self-interaction. */
if (gpart_is_ignorable(&gparts[i])) continue; fof_attach_self_cell(props, attach_r2, space_gparts, nr_gparts, c);
}
}
const size_t index = gparts[i].fof_data.group_id - group_id_offset; /* Mapper function to atomically update the group size array. */
void fof_update_group_size_mapper(hashmap_key_t key, hashmap_value_t *value,
void *data) {
/* Only check groups above the minimum mass threshold. */ size_t *group_size = (size_t *)data;
if (gparts[i].fof_data.group_id != group_id_default) {
/* Compute the centre of mass */ /* Use key to index into group size array. */
const double mass = gparts[i].mass; atomic_add(&group_size[key], value->value_st);
double x[3] = {gparts[i].x[0], gparts[i].x[1], gparts[i].x[2]}; }
/* Record the first particle of this group that we encounter so we /**
* can use it as reference frame for the centre of mass calculation */ * @brief Mapper function to calculate the group sizes.
if (first_position[index * 3 + 0] == (double)(-FLT_MAX)) { *
first_position[index * 3 + 0] = x[0]; * @param map_data An array of #gpart%s.
first_position[index * 3 + 1] = x[1]; * @param num_elements Chunk size.
first_position[index * 3 + 2] = x[2]; * @param extra_data Pointer to a #space.
} */
void fof_calc_group_size_mapper(void *map_data, int num_elements,
void *extra_data) {
if (periodic) { /* Retrieve mapped data. */
x[0] = nearest(x[0] - first_position[index * 3 + 0], dim[0]); struct space *s = (struct space *)extra_data;
x[1] = nearest(x[1] - first_position[index * 3 + 1], dim[1]); struct gpart *gparts = (struct gpart *)map_data;
x[2] = nearest(x[2] - first_position[index * 3 + 2], dim[2]); size_t *restrict group_index = s->e->fof_properties->group_index;
} size_t *restrict group_size = s->e->fof_properties->group_size;
centre_of_mass[index * 3 + 0] += mass * x[0]; /* Offset into gparts array. */
centre_of_mass[index * 3 + 1] += mass * x[1]; const ptrdiff_t gparts_offset = (ptrdiff_t)(gparts - s->gparts);
centre_of_mass[index * 3 + 2] += mass * x[2]; size_t *const group_index_offset = group_index + gparts_offset;
/* Check haloes above the seeding threshold */ /* Create hash table. */
if (group_mass[index] > seed_halo_mass) { hashmap_t map;
hashmap_init(&map);
/* Find the densest gas particle. for (int ind = 0; ind < num_elements; ind++) {
* Account for groups that already have a black hole and groups that
* contain no gas. */
if (gparts[i].type == swift_type_gas &&
max_part_density_index[index] != fof_halo_has_black_hole) {
const size_t gas_index = -gparts[i].id_or_neg_offset; const hashmap_key_t root =
const float rho_com = hydro_get_comoving_density(&parts[gas_index]); (hashmap_key_t)fof_find(group_index_offset[ind], group_index);
const size_t gpart_index = gparts_offset + ind;
/* Update index if a denser gas particle is found. */ /* Only add particles which aren't the root of a group. Stops groups of size
if (rho_com > max_part_density[index]) { * 1 being added to the hash table. */
max_part_density_index[index] = gas_index; if (root != gpart_index) {
max_part_density[index] = rho_com; hashmap_value_t *size = hashmap_get(&map, root);
}
}
/* If there is already a black hole in the group we don't need to
create a new one. */
else if (gparts[i].type == swift_type_black_hole) {
max_part_density_index[index] = fof_halo_has_black_hole;
}
} else { if (size != NULL)
max_part_density_index[index] = fof_halo_has_too_low_mass; (*size).value_st++;
} else
error("Couldn't find key (%zu) or create new one.", root);
} }
} }
props->extra_bh_seed_count = 0; /* Update the group size array. */
#endif if (map.size > 0)
hashmap_iterate(&map, fof_update_group_size_mapper, group_size);
hashmap_free(&map);
} }
/** /**
...@@ -2774,111 +2439,282 @@ void fof_find_foreign_links_mapper(void *map_data, int num_elements, ...@@ -2774,111 +2439,282 @@ void fof_find_foreign_links_mapper(void *map_data, int num_elements,
int *restrict group_link_count = &props->group_link_count; int *restrict group_link_count = &props->group_link_count;
struct fof_mpi **group_links = &props->group_links; struct fof_mpi **group_links = &props->group_links;
/* If the global group_links array is not big enough re-allocate it. */ /* If the global group_links array is not big enough re-allocate it. */
if (*group_link_count + local_link_count > *group_links_size) { if (*group_link_count + local_link_count > *group_links_size) {
const int old_size = *group_links_size;
const int new_size =
max(*group_link_count + local_link_count, 2 * old_size);
(*group_links) = (struct fof_mpi *)realloc(
*group_links, new_size * sizeof(struct fof_mpi));
*group_links_size = new_size;
message("Re-allocating global group links from %d to %d elements.",
old_size, new_size);
}
/* Copy the local links to the global list. */
for (int i = 0; i < local_link_count; i++) {
int found = 0;
/* Check that the links have not already been added to the list by another
* thread. */
for (int l = 0; l < *group_link_count; l++) {
if ((*group_links)[l].group_i == local_group_links[i].group_i &&
(*group_links)[l].group_j == local_group_links[i].group_j) {
found = 1;
break;
}
}
if (!found) {
(*group_links)[*group_link_count].group_i =
local_group_links[i].group_i;
(*group_links)[*group_link_count].group_i_size =
local_group_links[i].group_i_size;
(*group_links)[*group_link_count].group_j =
local_group_links[i].group_j;
(*group_links)[*group_link_count].group_j_size =
local_group_links[i].group_j_size;
(*group_link_count) = (*group_link_count) + 1;
}
}
}
/* Release lock. */
if (lock_unlock(&s->lock) != 0) error("Failed to unlock the space");
swift_free("fof_local_group_links", local_group_links);
#endif
}
/**
* @brief Compute the group properties for all groups.
*
* At the end of the process, all MPI ranks have all the information
* about every group. We compute:
* - Group size,
* - Group total mass,
* - Group centre of mass,
* - Maximal gas particle density,
* - Whether a group has a BH particle or not.
*
* TODO: Possible improvement is to delay the allocation of the CoM array
* until after the min/max density arrays have played their role.
* TODO: Can the loop over particles be threadpoolized?
*/
void fof_calc_group_mass(struct fof_props *props, const struct space *s,
int *restrict number_of_local_seeds,
int *restrict number_of_global_seeds) {
const size_t nr_gparts = s->nr_gparts;
const struct gpart *gparts = s->gparts;
const struct part *parts = s->parts;
const size_t group_id_default = props->group_id_default;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double seed_halo_mass = props->seed_halo_mass;
/* Direct pointers to the arrays */
long long *final_group_size = props->final_group_size;
double *group_mass = props->group_mass;
double *centre_of_mass = props->group_centre_of_mass;
char *has_black_hole = props->has_black_hole;
float *max_part_density = props->max_part_density;
/* Temporary arrays to help with the CoMs */
float *max_positions, *min_positions;
if (swift_memalign("fof_group_max_position", (void **)&max_positions, 32,
props->num_groups * 3 * sizeof(double)) != 0)
error("Unable to allocate memory for the max positions");
if (swift_memalign("fof_group_min_position", (void **)&min_positions, 32,
props->num_groups * 3 * sizeof(double)) != 0)
error("Unable to allocate memory for the min positions");
/* Initialise the min/max arrays to the limits */
for (size_t i = 0; i < 3 * (size_t)props->num_groups; ++i) {
min_positions[i] = DBL_MAX;
max_positions[i] = -DBL_MAX;
}
/* Collect information about the local particles and update the array of
* properties. Recall the array is as big as all the haloes accross
* all domains */
for (size_t i = 0; i < nr_gparts; i++) {
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
/* Check whether we ignore this particle type altogether */
if (gpart_is_ignorable(&gparts[i])) continue;
/* Ignore particles not in groups */
if (gparts[i].fof_data.group_id == group_id_default) continue;
if (gparts[i].fof_data.group_id > (size_t)props->num_groups)
error("Found an invalid group ID!");
/* Entry into the global list of group properties */
const size_t index = gparts[i].fof_data.group_id - 1;
/********************
* We know in which group this particle is: compute props
********************/
/* Count the number of particles */
final_group_size[index]++;
/* Add to the total mass */
group_mass[index] += gparts[i].mass;
/* Get the min/max position along each axis */
min_positions[index * 3 + 0] =
fmin(min_positions[index * 3 + 0], gparts[i].x[0]);
min_positions[index * 3 + 1] =
fmin(min_positions[index * 3 + 1], gparts[i].x[1]);
min_positions[index * 3 + 2] =
fmin(min_positions[index * 3 + 2], gparts[i].x[2]);
max_positions[index * 3 + 0] =
fmax(max_positions[index * 3 + 0], gparts[i].x[0]);
max_positions[index * 3 + 1] =
fmax(max_positions[index * 3 + 1], gparts[i].x[1]);
max_positions[index * 3 + 2] =
fmax(max_positions[index * 3 + 2], gparts[i].x[2]);
/* Check whether there is a black hole */
if (gparts[i].type == swift_type_black_hole) {
has_black_hole[index] = 1;
}
/* Idntify the densest gas particle in the group */
if (gparts[i].type == swift_type_gas) {
const size_t gas_index = -gparts[i].id_or_neg_offset;
const float rho_com = hydro_get_comoving_density(&parts[gas_index]);
max_part_density[index] = fmaxf(rho_com, max_part_density[index]);
}
}
const int old_size = *group_links_size; #ifdef WITH_MPI
const int new_size = /* Now all-reduce the local fragments so that everyone has a full catalog */
max(*group_link_count + local_link_count, 2 * old_size); MPI_Allreduce(MPI_IN_PLACE, final_group_size, props->num_groups,
MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(MPI_IN_PLACE, group_mass, props->num_groups, MPI_DOUBLE,
MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(MPI_IN_PLACE, min_positions, 3 * props->num_groups, MPI_DOUBLE,
MPI_MIN, MPI_COMM_WORLD);
MPI_Allreduce(MPI_IN_PLACE, max_positions, 3 * props->num_groups, MPI_DOUBLE,
MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(MPI_IN_PLACE, has_black_hole, props->num_groups, MPI_CHAR,
MPI_MAX, MPI_COMM_WORLD);
MPI_Allreduce(MPI_IN_PLACE, max_part_density, props->num_groups, MPI_FLOAT,
MPI_MAX, MPI_COMM_WORLD);
#endif
*number_of_local_seeds = 0;
*number_of_global_seeds = 0;
/* We can now do a second pass to compute the centre of mass
* Because of periodic BCs, we need to shift the positions in cases where the
* (max-min) is larger than 1/2 of the box size */
for (size_t i = 0; i < nr_gparts; i++) {
(*group_links) = (struct fof_mpi *)realloc( /* Ignore inhibited particles */
*group_links, new_size * sizeof(struct fof_mpi)); if (gparts[i].time_bin >= time_bin_inhibited) continue;
*group_links_size = new_size; /* Check whether we ignore this particle type altogether */
if (gpart_is_ignorable(&gparts[i])) continue;
message("Re-allocating global group links from %d to %d elements.", /* Ignore particles not in groups */
old_size, new_size); if (gparts[i].fof_data.group_id == group_id_default) continue;
}
/* Copy the local links to the global list. */ /* Entry into the global list of group properties */
for (int i = 0; i < local_link_count; i++) { const size_t index = gparts[i].fof_data.group_id - 1;
int found = 0; const int halo_is_on_edge[3] = {
max_positions[index * 3 + 0] - min_positions[index * 3 + 0] >
0.5 * dim[0],
max_positions[index * 3 + 1] - min_positions[index * 3 + 1] >
0.5 * dim[1],
max_positions[index * 3 + 2] - min_positions[index * 3 + 2] >
0.5 * dim[2]};
/* Check that the links have not already been added to the list by another /* Get particle position, including necessary wrapping */
* thread. */ double x[3] = {gparts[i].x[0], gparts[i].x[1], gparts[i].x[2]};
for (int l = 0; l < *group_link_count; l++) { if (periodic) {
if ((*group_links)[l].group_i == local_group_links[i].group_i && for (int k = 0; k < 3; k++) {
(*group_links)[l].group_j == local_group_links[i].group_j) { if (halo_is_on_edge[k]) {
found = 1; x[k] = box_wrap(x[k] + 0.5 * dim[k], 0., dim[k]);
break;
} }
} }
}
if (!found) { /* Centre of mass */
centre_of_mass[index * 3 + 0] += gparts[i].mass * x[0];
centre_of_mass[index * 3 + 1] += gparts[i].mass * x[1];
centre_of_mass[index * 3 + 2] += gparts[i].mass * x[2];
(*group_links)[*group_link_count].group_i = /* Should we seed a BH in this group? */
local_group_links[i].group_i; if (!has_black_hole[index] && group_mass[index] > seed_halo_mass) {
(*group_links)[*group_link_count].group_i_size =
local_group_links[i].group_i_size;
(*group_links)[*group_link_count].group_j = /* Is this a gas particle? */
local_group_links[i].group_j; if (gparts[i].type == swift_type_gas) {
(*group_links)[*group_link_count].group_j_size = const size_t gas_index = -gparts[i].id_or_neg_offset;
local_group_links[i].group_j_size; const float rho_com = hydro_get_comoving_density(&parts[gas_index]);
(*group_link_count) = (*group_link_count) + 1; /* Is this the gas paricle which is the densest? */
if (rho_com == max_part_density[index]) {
(*number_of_local_seeds)++;
}
} }
} }
} }
/* Release lock. */ #ifdef WITH_MPI
if (lock_unlock(&s->lock) != 0) error("Failed to unlock the space"); /* Now all-reduce the CoMs so that everyone has a full catalog */
MPI_Allreduce(MPI_IN_PLACE, centre_of_mass, 3 * props->num_groups, MPI_DOUBLE,
swift_free("fof_local_group_links", local_group_links); MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce(number_of_local_seeds, number_of_global_seeds, 1, MPI_INT,
MPI_SUM, MPI_COMM_WORLD);
#else
*number_of_global_seeds = *number_of_local_seeds;
#endif #endif
}
void fof_finalise_group_data(struct fof_props *props,
const struct group_length *group_sizes,
const struct gpart *gparts, const int periodic,
const double dim[3], const int num_groups) {
size_t *group_size = /* Finalise the operations on the group catalog */
(size_t *)swift_malloc("fof_group_size", num_groups * sizeof(size_t)); for (size_t i = 0; i < (size_t)props->num_groups; ++i) {
size_t *group_index =
(size_t *)swift_malloc("fof_group_index", num_groups * sizeof(size_t));
double *group_centre_of_mass = (double *)swift_malloc(
"fof_group_centre_of_mass", 3 * num_groups * sizeof(double));
for (int i = 0; i < num_groups; i++) { centre_of_mass[i * 3 + 0] /= group_mass[i];
centre_of_mass[i * 3 + 1] /= group_mass[i];
centre_of_mass[i * 3 + 2] /= group_mass[i];
const size_t group_offset = group_sizes[i].index; const int halo_is_on_edge[3] = {
max_positions[i * 3 + 0] - min_positions[i * 3 + 0] > 0.5 * dim[0],
max_positions[i * 3 + 1] - min_positions[i * 3 + 1] > 0.5 * dim[1],
max_positions[i * 3 + 2] - min_positions[i * 3 + 2] > 0.5 * dim[2]};
/* Centre of mass, including possible box wrapping */ /* Undo the half-box periodic shift */
double CoM[3] = {
props->group_centre_of_mass[i * 3 + 0] / props->group_mass[i],
props->group_centre_of_mass[i * 3 + 1] / props->group_mass[i],
props->group_centre_of_mass[i * 3 + 2] / props->group_mass[i]};
if (periodic) { if (periodic) {
CoM[0] = for (int k = 0; k < 3; k++) {
box_wrap(CoM[0] + props->group_first_position[i * 3 + 0], 0., dim[0]); if (halo_is_on_edge[k]) {
CoM[1] = centre_of_mass[i * 3 + k] -= 0.5 * dim[k];
box_wrap(CoM[1] + props->group_first_position[i * 3 + 1], 0., dim[1]); }
CoM[2] = centre_of_mass[i * 3 + k] =
box_wrap(CoM[2] + props->group_first_position[i * 3 + 2], 0., dim[2]); box_wrap(centre_of_mass[i * 3 + k], 0., dim[k]);
}
} }
#ifdef WITH_MPI /* Relabel the group ids */
group_index[i] = gparts[group_offset - node_offset].fof_data.group_id; props->final_group_index[i] = i + 1;
group_size[i] = props->group_size[group_offset - node_offset];
#else
group_index[i] = gparts[group_offset].fof_data.group_id;
group_size[i] = props->group_size[group_offset];
#endif
group_centre_of_mass[i * 3 + 0] = CoM[0];
group_centre_of_mass[i * 3 + 1] = CoM[1];
group_centre_of_mass[i * 3 + 2] = CoM[2];
} }
swift_free("fof_group_centre_of_mass", props->group_centre_of_mass); /* Free temporary arrays */
swift_free("fof_group_size", props->group_size); swift_free("max_positions", max_positions);
swift_free("fof_group_index", props->group_index); swift_free("min_positions", min_positions);
props->group_centre_of_mass = group_centre_of_mass;
props->group_size = group_size;
props->group_index = group_index;
} }
/** /**
...@@ -2890,41 +2726,25 @@ void fof_finalise_group_data(struct fof_props *props, ...@@ -2890,41 +2726,25 @@ void fof_finalise_group_data(struct fof_props *props,
* @param constants The physical constants. * @param constants The physical constants.
* @param cosmo The cosmological model. * @param cosmo The cosmological model.
* @param s The @space we act on. * @param s The @space we act on.
* @param num_groups_local The number of groups on the current MPI rank. * @param number_of_local_seeds Number of BHs to create on this rank.
* @param group_sizes List of groups sorted in size order. * @param number_of_global_seeds Number of BHs to create in total.
*/ */
void fof_seed_black_holes(const struct fof_props *props, void fof_seed_black_holes(const struct fof_props *props,
const struct black_holes_props *bh_props, const struct black_holes_props *bh_props,
const struct phys_const *constants, const struct phys_const *constants,
const struct cosmology *cosmo, struct space *s, const struct cosmology *cosmo, struct space *s,
const int num_groups_local) { const int number_of_local_seeds,
const int number_of_global_seeds) {
const long long *max_part_density_index = props->max_part_density_index;
/* Count the number of black holes to seed */
int num_seed_black_holes = 0;
for (int i = 0; i < num_groups_local + props->extra_bh_seed_count; i++) {
if (max_part_density_index[i] >= 0) ++num_seed_black_holes;
}
#ifdef WITH_MPI
int total_num_seed_black_holes = 0;
/* Sum the total number of black holes over each MPI rank. */
MPI_Reduce(&num_seed_black_holes, &total_num_seed_black_holes, 1, MPI_INT,
MPI_SUM, 0, MPI_COMM_WORLD);
#else
int total_num_seed_black_holes = num_seed_black_holes;
#endif
if (engine_rank == 0) if (engine_rank == 0)
message("Seeding %d black hole(s)", total_num_seed_black_holes); message("Seeding %d black hole(s)", number_of_global_seeds);
/* Anything to do this time on this rank? */ /* Anything to do this time on this rank? */
if (num_seed_black_holes == 0) return; if (number_of_local_seeds == 0) return;
/* Do we need to reallocate the black hole array for the new particles? */ /* Do we need to reallocate the black hole array for the new particles? */
if (s->nr_bparts + num_seed_black_holes > s->size_bparts) { if (s->nr_bparts + number_of_local_seeds > s->size_bparts) {
const size_t nr_bparts_new = s->nr_bparts + num_seed_black_holes; const size_t nr_bparts_new = s->nr_bparts + number_of_local_seeds;
s->size_bparts = engine_parts_size_grow * nr_bparts_new; s->size_bparts = engine_parts_size_grow * nr_bparts_new;
...@@ -2938,70 +2758,111 @@ void fof_seed_black_holes(const struct fof_props *props, ...@@ -2938,70 +2758,111 @@ void fof_seed_black_holes(const struct fof_props *props,
s->bparts = bparts_new; s->bparts = bparts_new;
} }
int k = s->nr_bparts; const size_t nr_gparts = s->nr_gparts;
struct gpart *gparts = s->gparts;
struct part *parts = s->parts;
struct xpart *xparts = s->xparts;
struct bpart *bparts = s->bparts;
const size_t group_id_default = props->group_id_default;
const double seed_halo_mass = props->seed_halo_mass;
/* Direct pointers to the arrays */
double *group_mass = props->group_mass;
char *has_black_hole = props->has_black_hole;
float *max_part_density = props->max_part_density;
size_t k = s->nr_bparts;
for (size_t i = 0; i < nr_gparts; i++) {
/* Check whether this is a gas particle */
if (gparts[i].type != swift_type_gas) continue;
/* Ignore inhibited particles */
if (gparts[i].time_bin >= time_bin_inhibited) continue;
/* Ignore particles not in groups */
if (gparts[i].fof_data.group_id == group_id_default) continue;
/* Loop over the local groups */ if (gparts[i].fof_data.group_id > (size_t)props->num_groups)
for (int i = 0; i < num_groups_local + props->extra_bh_seed_count; i++) { error("Found an invalid group ID!");
const long long part_index = max_part_density_index[i]; /* Get the density of this particle */
const size_t gas_index = -gparts[i].id_or_neg_offset;
const float rho_com = hydro_get_comoving_density(&parts[gas_index]);
/* Should we seed? */ /* Entry into the global list of group properties of this particle */
if (part_index >= 0) { const size_t index = gparts[i].fof_data.group_id - 1;
/* Handle on the particle to convert */ /* Should we seed a BH in this group? */
struct part *p = &s->parts[part_index]; if (!has_black_hole[index] && group_mass[index] > seed_halo_mass) {
struct xpart *xp = &s->xparts[part_index];
struct gpart *gp = p->gpart;
/* Let's destroy the gas particle */ /* Does it match the max density for this group?
p->time_bin = time_bin_inhibited; * (i.e. is it the particle we identified as the one to convert?) */
p->gpart = NULL; if (rho_com == max_part_density[index]) {
/* Mark the gpart as black hole */ /* Handle on the particle to convert */
gp->type = swift_type_black_hole; struct part *p = &parts[gas_index];
struct xpart *xp = &xparts[gas_index];
struct gpart *gp = p->gpart;
/* Basic properties of the black hole */ #ifdef SWIFT_DEBUG_CHECKS
struct bpart *bp = &s->bparts[k]; if (gp != &gparts[i]) error("Weird gas<->gpart linking error!");
bzero(bp, sizeof(struct bpart)); #endif
bp->time_bin = gp->time_bin;
/* Let's destroy the gas particle */
p->time_bin = time_bin_inhibited;
p->gpart = NULL;
/* Mark the gpart as black hole */
gp->type = swift_type_black_hole;
/* Re-link things */ /* Basic properties of the black hole */
bp->gpart = gp; struct bpart *bp = &bparts[k];
gp->id_or_neg_offset = -(bp - s->bparts); bzero(bp, sizeof(struct bpart));
bp->time_bin = gp->time_bin;
/* Synchronize masses, positions and velocities */ /* Re-link things */
bp->mass = gp->mass; bp->gpart = gp;
bp->x[0] = gp->x[0]; gp->id_or_neg_offset = -(bp - s->bparts);
bp->x[1] = gp->x[1];
bp->x[2] = gp->x[2];
bp->v[0] = gp->v_full[0];
bp->v[1] = gp->v_full[1];
bp->v[2] = gp->v_full[2];
/* Set a smoothing length */ /* Synchronize masses, positions and velocities */
bp->h = p->h; bp->mass = gp->mass;
bp->x[0] = gp->x[0];
bp->x[1] = gp->x[1];
bp->x[2] = gp->x[2];
bp->v[0] = gp->v_full[0];
bp->v[1] = gp->v_full[1];
bp->v[2] = gp->v_full[2];
/* Save the ID */ /* Set a smoothing length */
bp->id = p->id; bp->h = p->h;
/* Save the ID */
bp->id = p->id;
/* Save the tree depth */
bp->depth_h = p->depth_h;
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
bp->ti_kick = p->ti_kick; bp->ti_kick = p->ti_kick;
bp->ti_drift = p->ti_drift; bp->ti_drift = p->ti_drift;
#endif #endif
/* Copy over all the gas properties that we want */ /* Copy over all the gas properties that we want */
black_holes_create_from_gas(bp, bh_props, constants, cosmo, p, xp, black_holes_create_from_gas(bp, bh_props, constants, cosmo, p, xp,
s->e->ti_current); s->e->ti_current);
tracers_first_init_bpart(bp, s->e->internal_units, tracers_first_init_bpart(bp, s->e->internal_units,
s->e->physical_constants, cosmo); s->e->physical_constants, cosmo);
/* Move to the next BH slot */ /* Move to the next BH slot */
k++; k++;
}
} }
} }
#ifdef SWIFT_DEBUG_CHECKS #ifdef SWIFT_DEBUG_CHECKS
if ((int)(s->nr_bparts) + num_seed_black_holes != k) { if (s->nr_bparts + number_of_local_seeds != k) {
error("Seeded the wrong number of black holes!"); error("Seeded the wrong number of black holes!");
} }
#endif #endif
...@@ -3017,21 +2878,14 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank, ...@@ -3017,21 +2878,14 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank,
FILE *file = NULL; FILE *file = NULL;
struct part *parts = s->parts;
long long *final_group_size = props->final_group_size; long long *final_group_size = props->final_group_size;
size_t *group_index = props->group_index; long long *group_index = props->final_group_index;
double *group_mass = props->group_mass; double *group_mass = props->group_mass;
double *group_centre_of_mass = props->group_centre_of_mass; double *group_centre_of_mass = props->group_centre_of_mass;
const long long *max_part_density_index = props->max_part_density_index;
const float *max_part_density = props->max_part_density;
for (int rank = 0; rank < nr_nodes; ++rank) { for (int rank = 0; rank < nr_nodes; ++rank) {
#ifdef WITH_MPI if (rank == 0) {
MPI_Barrier(MPI_COMM_WORLD);
#endif
if (rank == my_rank) {
const char *mode; const char *mode;
if (my_rank == 0) if (my_rank == 0)
...@@ -3058,27 +2912,12 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank, ...@@ -3058,27 +2912,12 @@ void fof_dump_group_data(const struct fof_props *props, const int my_rank,
for (int i = 0; i < num_groups; i++) { for (int i = 0; i < num_groups; i++) {
const long long part_id = props->max_part_density_index[i] >= 0 fprintf(file, " %8lld %12lld %12e %12e %12e %12e %12e %24lld %24lld\n",
? parts[max_part_density_index[i]].id
: -1;
fprintf(file, " %8zu %12lld %12e %12e %12e %12e %12e %24lld %24lld\n",
group_index[i], final_group_size[i], group_mass[i], group_index[i], final_group_size[i], group_mass[i],
group_centre_of_mass[i * 3 + 0], group_centre_of_mass[i * 3 + 0],
group_centre_of_mass[i * 3 + 1], group_centre_of_mass[i * 3 + 1],
group_centre_of_mass[i * 3 + 2], max_part_density[i], group_centre_of_mass[i * 3 + 2], 0., -1ll, -1ll);
max_part_density_index[i], part_id);
}
/* Dump the extra black hole seeds. */
for (int i = num_groups; i < num_groups + props->extra_bh_seed_count;
i++) {
const long long part_id = max_part_density_index[i] >= 0
? parts[max_part_density_index[i]].id
: -1;
fprintf(file, " %8zu %12zu %12e %12e %12e %12e %12e %24lld %24lld\n",
0UL, 0UL, 0., 0., 0., 0., 0., 0LL, part_id);
} }
fclose(file); fclose(file);
} }
} }
...@@ -3396,90 +3235,6 @@ void fof_link_attachable_particles(struct fof_props *props, ...@@ -3396,90 +3235,6 @@ void fof_link_attachable_particles(struct fof_props *props,
clocks_from_ticks(getticks() - tic_total), clocks_getunit()); clocks_from_ticks(getticks() - tic_total), clocks_getunit());
} }
/**
* @brief Process all the attachable-linkable connections to add the
* attachables to the groups they belong to.
*
* @param props The properties fof the FOF scheme.
* @param s The #space we work with.
*/
void fof_finalise_attachables(struct fof_props *props, const struct space *s) {
/* Is there anything to attach? */
if (!current_fof_attach_type) return;
const ticks tic_total = getticks();
const size_t nr_gparts = s->nr_gparts;
char *found_attachable_link = props->found_attachable_link;
size_t *restrict group_index = props->group_index;
size_t *restrict attach_index = props->attach_index;
size_t *restrict group_size = props->group_size;
#ifdef WITH_MPI
/* Get pointers to global arrays. */
int *restrict group_links_size = &props->group_links_size;
int *restrict group_link_count = &props->group_link_count;
struct fof_mpi **group_links = &props->group_links;
/* Loop over all the attachables and added them to the group they belong to */
for (size_t i = 0; i < nr_gparts; ++i) {
const struct gpart *gp = &s->gparts[i];
if (gpart_is_attachable(gp) && found_attachable_link[i]) {
/* Update its root */
const size_t root_j = attach_index[i];
const size_t root_i =
fof_find_global(group_index[i] - node_offset, group_index, nr_gparts);
/* Update the size of the group the particle belongs to */
if (is_local(root_j, nr_gparts)) {
const size_t local_root = root_j - node_offset;
group_index[i] = local_root + node_offset;
group_size[local_root]++;
} else {
add_foreign_link_to_list(group_link_count, group_links_size,
group_links, group_links, root_i, root_j,
/*size_i=*/1,
/*size_j=*/2);
}
}
}
#else /* not WITH_MPI */
/* Loop over all the attachables and added them to the group they belong to */
for (size_t i = 0; i < nr_gparts; ++i) {
const struct gpart *gp = &s->gparts[i];
if (gpart_is_attachable(gp) && found_attachable_link[i]) {
const size_t root = attach_index[i];
/* Update its root */
group_index[i] = root;
/* Update the size of the group the particle belongs to */
group_size[root]++;
}
}
#endif /* WITH_MPI */
if (s->e->verbose)
message("fof_finalise_attachables() took (FOF SCALING): %.3f %s.",
clocks_from_ticks(getticks() - tic_total), clocks_getunit());
}
/** /**
* @brief Process all the group fragments spanning more than * @brief Process all the group fragments spanning more than
* one rank to link them. * one rank to link them.
...@@ -3522,6 +3277,9 @@ void fof_link_foreign_fragments(struct fof_props *props, ...@@ -3522,6 +3277,9 @@ void fof_link_foreign_fragments(struct fof_props *props,
MPI_Allreduce(&group_link_count, &global_group_link_count, 1, MPI_INT, MPI_Allreduce(&group_link_count, &global_group_link_count, 1, MPI_INT,
MPI_SUM, MPI_COMM_WORLD); MPI_SUM, MPI_COMM_WORLD);
if (global_group_link_count < 0)
error("Overflow of the size of the global list of foreign links");
struct fof_mpi *global_group_links = NULL; struct fof_mpi *global_group_links = NULL;
int *displ = NULL, *group_link_counts = NULL; int *displ = NULL, *group_link_counts = NULL;
...@@ -3807,13 +3565,7 @@ void fof_compute_local_sizes(struct fof_props *props, struct space *s) { ...@@ -3807,13 +3565,7 @@ void fof_compute_local_sizes(struct fof_props *props, struct space *s) {
* @param dump_results Do we want to write the group catalogue to a hdf5 file? * @param dump_results Do we want to write the group catalogue to a hdf5 file?
* @param seed_black_holes Do we want to seed black holes in haloes? * @param seed_black_holes Do we want to seed black holes in haloes?
*/ */
void fof_compute_group_props(struct fof_props *props, void fof_assign_group_ids(struct fof_props *props, struct space *s) {
const struct black_holes_props *bh_props,
const struct phys_const *constants,
const struct cosmology *cosmo, struct space *s,
const int dump_results,
const int dump_debug_results,
const int seed_black_holes) {
const int verbose = s->e->verbose; const int verbose = s->e->verbose;
#ifdef WITH_MPI #ifdef WITH_MPI
...@@ -3868,7 +3620,7 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3868,7 +3620,7 @@ void fof_compute_group_props(struct fof_props *props,
/* Sort the groups in descending order based upon size and re-label their /* Sort the groups in descending order based upon size and re-label their
* IDs 0-num_groups. */ * IDs 0-num_groups. */
struct group_length *high_group_sizes = NULL; struct group_length *high_group_sizes = NULL;
int group_count = 0; size_t group_count = 0;
if (swift_memalign("fof_high_group_sizes", (void **)&high_group_sizes, 32, if (swift_memalign("fof_high_group_sizes", (void **)&high_group_sizes, 32,
num_groups_local * sizeof(struct group_length)) != 0) num_groups_local * sizeof(struct group_length)) != 0)
...@@ -3889,6 +3641,7 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3889,6 +3641,7 @@ void fof_compute_group_props(struct fof_props *props,
} }
#endif #endif
} }
props->high_group_sizes = high_group_sizes;
ticks tic = getticks(); ticks tic = getticks();
...@@ -3897,30 +3650,25 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3897,30 +3650,25 @@ void fof_compute_group_props(struct fof_props *props,
#ifdef WITH_MPI #ifdef WITH_MPI
MPI_Allreduce(&num_groups_local, &num_groups, 1, MPI_LONG_LONG_INT, MPI_SUM, MPI_Allreduce(&num_groups_local, &num_groups, 1, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD); MPI_COMM_WORLD);
if (verbose)
message("Finding the total no. of groups took: (FOF SCALING): %.3f %s.",
clocks_from_ticks(getticks() - tic_num_groups_calc),
clocks_getunit());
MPI_Reduce(&num_parts_in_groups_local, &num_parts_in_groups, 1, MPI_Reduce(&num_parts_in_groups_local, &num_parts_in_groups, 1,
MPI_LONG_LONG_INT, MPI_SUM, 0, MPI_COMM_WORLD); MPI_LONG_LONG_INT, MPI_SUM, 0, MPI_COMM_WORLD);
MPI_Reduce(&max_group_size_local, &max_group_size, 1, MPI_LONG_LONG_INT, MPI_Reduce(&max_group_size_local, &max_group_size, 1, MPI_LONG_LONG_INT,
MPI_MAX, 0, MPI_COMM_WORLD); MPI_MAX, 0, MPI_COMM_WORLD);
#else #else
num_groups = num_groups_local; num_groups = num_groups_local;
num_parts_in_groups = num_parts_in_groups_local; num_parts_in_groups = num_parts_in_groups_local;
max_group_size = max_group_size_local; max_group_size = max_group_size_local;
#endif /* WITH_MPI */ #endif /* WITH_MPI */
props->num_groups = num_groups; props->num_groups = num_groups;
// message("num_groups_local=%zd", num_groups_local);
/* Find number of groups on lower numbered MPI ranks */ /* Find number of groups on lower numbered MPI ranks */
#ifdef WITH_MPI #ifdef WITH_MPI
long long nglocal = num_groups_local; long long nglocal = num_groups_local;
long long ngsum; long long ngsum;
MPI_Scan(&nglocal, &ngsum, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD); MPI_Scan(&nglocal, &ngsum, 1, MPI_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
const size_t num_groups_prev = (size_t)(ngsum - nglocal); props->num_groups_prev = (size_t)(ngsum - nglocal);
#endif /* WITH_MPI */ #endif /* WITH_MPI */
if (verbose) if (verbose)
...@@ -3949,7 +3697,7 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3949,7 +3697,7 @@ void fof_compute_group_props(struct fof_props *props,
for (size_t i = 0; i < num_groups_local; i++) { for (size_t i = 0; i < num_groups_local; i++) {
#ifdef WITH_MPI #ifdef WITH_MPI
gparts[high_group_sizes[i].index - node_offset].fof_data.group_id = gparts[high_group_sizes[i].index - node_offset].fof_data.group_id =
group_id_offset + i + num_groups_prev; group_id_offset + i + props->num_groups_prev;
#else #else
gparts[high_group_sizes[i].index].fof_data.group_id = group_id_offset + i; gparts[high_group_sizes[i].index].fof_data.group_id = group_id_offset + i;
#endif #endif
...@@ -3961,7 +3709,7 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3961,7 +3709,7 @@ void fof_compute_group_props(struct fof_props *props,
* AND the total size of the group is >= min_group_size we need to * AND the total size of the group is >= min_group_size we need to
* retrieve the gparts.group_id we just assigned to the global root. * retrieve the gparts.group_id we just assigned to the global root.
* *
* Will do that by sending the group_index of these lcoal roots to the * Will do that by sending the group_index of these local roots to the
* node where their global root is stored and receiving back the new * node where their global root is stored and receiving back the new
* group_id associated with that particle. * group_id associated with that particle.
* *
...@@ -3997,22 +3745,24 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -3997,22 +3745,24 @@ void fof_compute_group_props(struct fof_props *props,
compare_fof_final_index_global_root); compare_fof_final_index_global_root);
/* Determine range of global indexes (i.e. particles) on each node */ /* Determine range of global indexes (i.e. particles) on each node */
size_t *num_on_node = (size_t *)malloc(nr_nodes * sizeof(size_t)); props->num_on_node =
MPI_Allgather(&nr_gparts, sizeof(size_t), MPI_BYTE, num_on_node, (size_t *)swift_malloc("fof_num_on_node", nr_nodes * sizeof(size_t));
MPI_Allgather(&nr_gparts, sizeof(size_t), MPI_BYTE, props->num_on_node,
sizeof(size_t), MPI_BYTE, MPI_COMM_WORLD); sizeof(size_t), MPI_BYTE, MPI_COMM_WORLD);
size_t *first_on_node = (size_t *)malloc(nr_nodes * sizeof(size_t)); props->first_on_node =
first_on_node[0] = 0; (size_t *)swift_malloc("fof_first_on_node", nr_nodes * sizeof(size_t));
for (int i = 1; i < nr_nodes; i += 1) props->first_on_node[0] = 0;
first_on_node[i] = first_on_node[i - 1] + num_on_node[i - 1]; for (int i = 1; i < nr_nodes; i++)
props->first_on_node[i] =
props->first_on_node[i - 1] + props->num_on_node[i - 1];
/* Determine how many entries go to each node */ /* Determine how many entries go to each node */
int *sendcount = (int *)malloc(nr_nodes * sizeof(int)); int *sendcount = (int *)calloc(nr_nodes, sizeof(int));
for (int i = 0; i < nr_nodes; i += 1) sendcount[i] = 0;
int dest = 0; int dest = 0;
for (size_t i = 0; i < nsend; i += 1) { for (size_t i = 0; i < nsend; i += 1) {
while ((fof_index_send[i].global_root >= while ((fof_index_send[i].global_root >=
first_on_node[dest] + num_on_node[dest]) || props->first_on_node[dest] + props->num_on_node[dest]) ||
(num_on_node[dest] == 0)) (props->num_on_node[dest] == 0))
dest += 1; dest += 1;
if (dest >= nr_nodes) error("Node index out of range!"); if (dest >= nr_nodes) error("Node index out of range!");
sendcount[dest] += 1; sendcount[dest] += 1;
...@@ -4068,91 +3818,119 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -4068,91 +3818,119 @@ void fof_compute_group_props(struct fof_props *props,
#endif /* WITH_MPI */ #endif /* WITH_MPI */
size_t max_id = 0;
/* Assign every particle the group_id of its local root. */ /* Assign every particle the group_id of its local root. */
for (size_t i = 0; i < nr_gparts; i++) { for (size_t i = 0; i < nr_gparts; i++) {
if (gpart_is_ignorable(&gparts[i])) continue;
if (gpart_is_attachable(&gparts[i])) continue;
const size_t root = fof_find_local(i, nr_gparts, group_index); const size_t root = fof_find_local(i, nr_gparts, group_index);
gparts[i].fof_data.group_id = gparts[root].fof_data.group_id; gparts[i].fof_data.group_id = gparts[root].fof_data.group_id;
if (gparts[i].fof_data.group_id != fof_props_default_group_id)
max_id = max(max_id, gparts[i].fof_data.group_id);
}
/* Give some info */
if (engine_rank == 0) {
message(
"No. of groups: %lld. No. of particles in groups: %lld. No. of "
"particles not in groups: %lld.",
num_groups, num_parts_in_groups,
s->e->total_nr_gparts - num_parts_in_groups);
message("Largest group (linkables only) by size: %lld", max_group_size);
} }
/* Free data we are done with */
swift_free("fof_group_index", props->group_index);
props->group_index = NULL;
if (verbose) if (verbose)
message("Group sorting took: %.3f %s.", clocks_from_ticks(getticks() - tic), message("took: %.3f %s.", clocks_from_ticks(getticks() - tic_total),
clocks_getunit()); clocks_getunit());
}
/**
* @brief Compute all the group properties
*
* @param props The properties of the FOF scheme.
* @param bh_props The properties of the black hole scheme.
* @param constants The physical constants in internal units.
* @param cosmo The current cosmological model.
* @param s The #space containing the particles.
* @param dump_debug_results Are we writing txt-file debug catalogues including
* BH-seeding info?
* @param dump_results Do we want to write the group catalogue to a hdf5 file?
* @param seed_black_holes Do we want to seed black holes in haloes?
*/
void fof_compute_group_props(struct fof_props *props,
const struct black_holes_props *bh_props,
const struct phys_const *constants,
const struct cosmology *cosmo, struct space *s,
const int dump_results,
const int dump_debug_results,
const int seed_black_holes) {
const int verbose = s->e->verbose;
const ticks tic_total = getticks();
/* Allocate and initialise a group mass and centre of mass array. */ const size_t num_groups = props->num_groups;
/* Allocate and initialise a group mass and centre of mass array
for *all* groups. */
if (swift_memalign("fof_group_mass", (void **)&props->group_mass, 32, if (swift_memalign("fof_group_mass", (void **)&props->group_mass, 32,
num_groups_local * sizeof(double)) != 0) num_groups * sizeof(double)) != 0)
error("Failed to allocate list of group masses for FOF search."); error("Failed to allocate list of group masses for FOF search.");
if (swift_memalign("fof_group_size", (void **)&props->final_group_size, 32, if (swift_memalign("fof_group_size", (void **)&props->final_group_size, 32,
num_groups_local * sizeof(long long)) != 0) num_groups * sizeof(long long)) != 0)
error("Failed to allocate list of group masses for FOF search.");
if (swift_memalign("fof_group_index", (void **)&props->final_group_index, 32,
num_groups * sizeof(long long)) != 0)
error("Failed to allocate list of group masses for FOF search."); error("Failed to allocate list of group masses for FOF search.");
if (swift_memalign("fof_has_black_hole", (void **)&props->has_black_hole, 32,
num_groups * sizeof(char)) != 0)
error("Failed to allocate list of black holes for FOF search.");
if (swift_memalign("fof_group_centre_of_mass", if (swift_memalign("fof_group_centre_of_mass",
(void **)&props->group_centre_of_mass, 32, (void **)&props->group_centre_of_mass, 32,
num_groups_local * 3 * sizeof(double)) != 0) num_groups * 3 * sizeof(double)) != 0)
error("Failed to allocate list of group CoM for FOF search."); error("Failed to allocate list of group CoM for FOF search.");
if (swift_memalign("fof_group_first_position",
(void **)&props->group_first_position, 32,
num_groups_local * 3 * sizeof(double)) != 0)
error("Failed to allocate list of group first positions for FOF search.");
bzero(props->group_mass, num_groups_local * sizeof(double));
bzero(props->final_group_size, num_groups_local * sizeof(long long));
bzero(props->group_centre_of_mass, num_groups_local * 3 * sizeof(double));
for (size_t i = 0; i < 3 * num_groups_local; i++) {
props->group_first_position[i] = -FLT_MAX;
}
/* Allocate and initialise arrays to identify the densest gas particle. */
if (swift_memalign("fof_max_part_density_index",
(void **)&props->max_part_density_index, 32,
num_groups_local * sizeof(long long)) != 0)
error(
"Failed to allocate list of max group density indices for FOF "
"search.");
if (swift_memalign("fof_max_part_density", (void **)&props->max_part_density, if (swift_memalign("fof_max_part_density", (void **)&props->max_part_density,
32, num_groups_local * sizeof(float)) != 0) 32, num_groups * sizeof(float)) != 0)
error("Failed to allocate list of max group densities for FOF search."); error("Failed to allocate list of max group densities for FOF search.");
/* No densest particle found so far */ bzero(props->group_mass, num_groups * sizeof(double));
bzero(props->max_part_density, num_groups_local * sizeof(float)); bzero(props->final_group_size, num_groups * sizeof(long long));
bzero(props->final_group_index, num_groups * sizeof(long long));
bzero(props->has_black_hole, num_groups * sizeof(char));
bzero(props->group_centre_of_mass, num_groups * 3 * sizeof(double));
bzero(props->max_part_density, num_groups * sizeof(float));
/* Start by assuming that the haloes have no gas */ const ticks tic_props = getticks();
for (size_t i = 0; i < num_groups_local; i++) {
props->max_part_density_index[i] = fof_halo_has_no_gas;
}
const ticks tic_seeding = getticks();
#ifdef WITH_MPI
fof_calc_group_mass(props, s, seed_black_holes, num_groups_local,
num_groups_prev, num_on_node, first_on_node,
props->group_mass);
free(num_on_node);
free(first_on_node);
#else
fof_calc_group_mass(props, s, seed_black_holes, num_groups_local,
/*num_groups_prev=*/0, /*num_on_node=*/NULL,
/*first_on_node=*/NULL, props->group_mass);
#endif
/* Finalise the group data before dump */ /* Now, compute all things */
fof_finalise_group_data(props, high_group_sizes, s->gparts, s->periodic, int number_of_local_seeds = 0, number_of_global_seeds = 0;
s->dim, num_groups_local); fof_calc_group_mass(props, s, &number_of_local_seeds,
&number_of_global_seeds);
if (verbose) if (verbose)
message("Computing group properties took: %.3f %s.", message("Computing group properties took: %.3f %s.",
clocks_from_ticks(getticks() - tic_seeding), clocks_getunit()); clocks_from_ticks(getticks() - tic_props), clocks_getunit());
/* All MPI ranks now have information about all haloes, even
* the ones where there are no local particles in. */
/* Dump group data. */ /* Dump group data (only rank 0 since everyone has everything anyway). */
if (dump_results) { if (dump_results && engine_rank == 0) {
#ifdef HAVE_HDF5 #ifdef HAVE_HDF5
write_fof_hdf5_catalogue(props, num_groups_local, s->e); write_fof_hdf5_catalogue(props, s->e);
#else #else
error("Can't dump hdf5 catalogues with hdf5 switched off!"); error("Can't dump hdf5 catalogues with hdf5 switched off!");
#endif #endif
} }
if (dump_debug_results) { if (dump_debug_results && engine_rank == 0) {
char output_file_name[PARSER_MAX_LINE_SIZE]; char output_file_name[PARSER_MAX_LINE_SIZE];
snprintf(output_file_name, PARSER_MAX_LINE_SIZE, "%s", props->base_name); snprintf(output_file_name, PARSER_MAX_LINE_SIZE, "%s", props->base_name);
...@@ -4165,52 +3943,49 @@ void fof_compute_group_props(struct fof_props *props, ...@@ -4165,52 +3943,49 @@ void fof_compute_group_props(struct fof_props *props,
".dat"); ".dat");
#endif #endif
fof_dump_group_data(props, s->e->nodeID, s->e->nr_nodes, output_file_name, fof_dump_group_data(props, s->e->nodeID, s->e->nr_nodes, output_file_name,
s, num_groups_local); s, num_groups);
} }
/* Seed black holes */ /* Seed black holes */
if (seed_black_holes) { if (seed_black_holes) {
fof_seed_black_holes(props, bh_props, constants, cosmo, s, fof_seed_black_holes(props, bh_props, constants, cosmo, s,
num_groups_local); number_of_local_seeds, number_of_global_seeds);
} }
/* Free the left-overs */ if (verbose)
swift_free("fof_high_group_sizes", high_group_sizes); message("took %.3f %s.", clocks_from_ticks(getticks() - tic_total),
clocks_getunit());
}
/**
* @brief Free all the arrays we allocated on the way
*/
void fof_free_arrays(struct fof_props *props) {
swift_free("fof_high_group_sizes", props->high_group_sizes);
swift_free("fof_group_mass", props->group_mass); swift_free("fof_group_mass", props->group_mass);
swift_free("fof_group_size", props->final_group_size); swift_free("fof_group_size", props->final_group_size);
swift_free("fof_group_index", props->final_group_index);
swift_free("fof_group_centre_of_mass", props->group_centre_of_mass); swift_free("fof_group_centre_of_mass", props->group_centre_of_mass);
swift_free("fof_group_first_position", props->group_first_position);
swift_free("fof_max_part_density_index", props->max_part_density_index);
swift_free("fof_max_part_density", props->max_part_density); swift_free("fof_max_part_density", props->max_part_density);
props->group_mass = NULL; swift_free("fof_has_black_hole", props->has_black_hole);
props->final_group_size = NULL;
props->group_centre_of_mass = NULL;
props->max_part_density_index = NULL;
props->max_part_density = NULL;
swift_free("fof_distance", props->distance_to_link); swift_free("fof_distance", props->distance_to_link);
swift_free("fof_group_index", props->group_index);
swift_free("fof_attach_index", props->attach_index); swift_free("fof_attach_index", props->attach_index);
swift_free("fof_found_attach", props->found_attachable_link); swift_free("fof_found_attach", props->found_attachable_link);
swift_free("fof_group_size", props->group_size); swift_free("fof_group_size", props->group_size);
props->group_index = NULL; props->group_mass = NULL;
props->final_group_size = NULL;
props->final_group_index = NULL;
props->group_centre_of_mass = NULL;
props->max_part_density = NULL;
props->has_black_hole = NULL;
props->group_size = NULL; props->group_size = NULL;
if (engine_rank == 0) {
message(
"No. of groups: %lld. No. of particles in groups: %lld. No. of "
"particles not in groups: %lld.",
num_groups, num_parts_in_groups,
s->e->total_nr_gparts - num_parts_in_groups);
message("Largest group by size: %lld", max_group_size);
}
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic_total),
clocks_getunit());
#ifdef WITH_MPI #ifdef WITH_MPI
MPI_Barrier(MPI_COMM_WORLD); swift_free("fof_num_on_node", props->num_on_node);
swift_free("fof_first_on_node", props->first_on_node);
props->num_on_node = NULL;
props->first_on_node = NULL;
#endif #endif
} }
...@@ -4225,7 +4000,6 @@ void fof_struct_dump(const struct fof_props *props, FILE *stream) { ...@@ -4225,7 +4000,6 @@ void fof_struct_dump(const struct fof_props *props, FILE *stream) {
temp.group_mass = NULL; temp.group_mass = NULL;
temp.final_group_size = NULL; temp.final_group_size = NULL;
temp.group_centre_of_mass = NULL; temp.group_centre_of_mass = NULL;
temp.max_part_density_index = NULL;
temp.max_part_density = NULL; temp.max_part_density = NULL;
temp.group_links = NULL; temp.group_links = NULL;
......
...@@ -37,6 +37,13 @@ struct phys_const; ...@@ -37,6 +37,13 @@ struct phys_const;
struct black_holes_props; struct black_holes_props;
struct cosmology; struct cosmology;
/* Store group size and offset into array. */
struct group_length {
size_t index, size;
} SWIFT_STRUCT_ALIGN;
struct fof_props { struct fof_props {
/*! Whether we're doing periodic FoF calls to seed black holes. */ /*! Whether we're doing periodic FoF calls to seed black holes. */
...@@ -80,10 +87,6 @@ struct fof_props { ...@@ -80,10 +87,6 @@ struct fof_props {
/*! Number of groups */ /*! Number of groups */
long long num_groups; long long num_groups;
/*! Number of local black holes that belong to groups whose roots are on a
* different node. */
int extra_bh_seed_count;
/*! Index of the root particle of the group a given gpart belongs to. */ /*! Index of the root particle of the group a given gpart belongs to. */
size_t *group_index; size_t *group_index;
...@@ -99,24 +102,36 @@ struct fof_props { ...@@ -99,24 +102,36 @@ struct fof_props {
/*! Size of the group a given gpart belongs to. */ /*! Size of the group a given gpart belongs to. */
size_t *group_size; size_t *group_size;
/*! Size of the local groups a given gpart belongs to. */
struct group_length *high_group_sizes;
/*! Final size of the group a given gpart belongs to. */ /*! Final size of the group a given gpart belongs to. */
long long *final_group_size; long long *final_group_size;
/*! Final index of the group a given gpart belongs to. */
long long *final_group_index;
/*! Mass of the group a given gpart belongs to. */ /*! Mass of the group a given gpart belongs to. */
double *group_mass; double *group_mass;
/*! Does the group have a black hole? */
char *has_black_hole;
/*! Centre of mass of the group a given gpart belongs to. */ /*! Centre of mass of the group a given gpart belongs to. */
double *group_centre_of_mass; double *group_centre_of_mass;
/*! Position of the first particle of a given group. */
double *group_first_position;
/*! Index of the part with the maximal density of each group. */
long long *max_part_density_index;
/*! Maximal density of all parts of each group. */ /*! Maximal density of all parts of each group. */
float *max_part_density; float *max_part_density;
/* Number of groups on each node */
size_t *num_on_node;
/* First group on each node */
size_t *first_on_node;
/* Total number of groups on lower numbered MPI ranks */
size_t num_groups_prev;
/* ------------ MPI-related arrays --------------- */ /* ------------ MPI-related arrays --------------- */
/*! The number of links between pairs of particles on this node and /*! The number of links between pairs of particles on this node and
...@@ -131,13 +146,6 @@ struct fof_props { ...@@ -131,13 +146,6 @@ struct fof_props {
struct fof_mpi *group_links; struct fof_mpi *group_links;
}; };
/* Store group size and offset into array. */
struct group_length {
size_t index, size;
} SWIFT_STRUCT_ALIGN;
#ifdef WITH_MPI #ifdef WITH_MPI
/* MPI message required for FOF. */ /* MPI message required for FOF. */
...@@ -196,8 +204,9 @@ void fof_compute_local_sizes(struct fof_props *props, struct space *s); ...@@ -196,8 +204,9 @@ void fof_compute_local_sizes(struct fof_props *props, struct space *s);
void fof_search_foreign_cells(struct fof_props *props, const struct space *s); void fof_search_foreign_cells(struct fof_props *props, const struct space *s);
void fof_link_attachable_particles(struct fof_props *props, void fof_link_attachable_particles(struct fof_props *props,
const struct space *s); const struct space *s);
void fof_finalise_attachables(struct fof_props *props, const struct space *s); void fof_finalise_attachables(struct fof_props *props, struct space *s);
void fof_link_foreign_fragments(struct fof_props *props, const struct space *s); void fof_link_foreign_fragments(struct fof_props *props, const struct space *s);
void fof_assign_group_ids(struct fof_props *props, struct space *s);
void fof_compute_group_props(struct fof_props *props, void fof_compute_group_props(struct fof_props *props,
const struct black_holes_props *bh_props, const struct black_holes_props *bh_props,
const struct phys_const *constants, const struct phys_const *constants,
...@@ -223,6 +232,7 @@ void rec_fof_attach_pair(const struct fof_props *props, const double dim[3], ...@@ -223,6 +232,7 @@ void rec_fof_attach_pair(const struct fof_props *props, const double dim[3],
const size_t nr_gparts, struct cell *restrict ci, const size_t nr_gparts, struct cell *restrict ci,
struct cell *restrict cj, const int ci_local, struct cell *restrict cj, const int ci_local,
const int cj_local); const int cj_local);
void fof_free_arrays(struct fof_props *props);
void fof_struct_dump(const struct fof_props *props, FILE *stream); void fof_struct_dump(const struct fof_props *props, FILE *stream);
void fof_struct_restore(struct fof_props *props, FILE *stream); void fof_struct_restore(struct fof_props *props, FILE *stream);
......
...@@ -64,9 +64,6 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e, ...@@ -64,9 +64,6 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
/* We write rank 0's hostname so that it is uniform across all files. */ /* We write rank 0's hostname so that it is uniform across all files. */
char systemname[256] = {0}; char systemname[256] = {0};
if (e->nodeID == 0) sprintf(systemname, "%s", hostname()); if (e->nodeID == 0) sprintf(systemname, "%s", hostname());
#ifdef WITH_MPI
if (!virtual_file) MPI_Bcast(systemname, 256, MPI_CHAR, 0, MPI_COMM_WORLD);
#endif
io_write_attribute_s(h_grp, "System", systemname); io_write_attribute_s(h_grp, "System", systemname);
io_write_attribute(h_grp, "Shift", DOUBLE, e->s->initial_shift, 3); io_write_attribute(h_grp, "Shift", DOUBLE, e->s->initial_shift, 3);
...@@ -115,8 +112,8 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e, ...@@ -115,8 +112,8 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
flagEntropy[0] = writeEntropyFlag(); flagEntropy[0] = writeEntropyFlag();
io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy, io_write_attribute(h_grp, "Flag_Entropy_ICs", UINT, flagEntropy,
swift_type_count); swift_type_count);
io_write_attribute_i(h_grp, "NumFilesPerSnapshot", e->nr_nodes); io_write_attribute_i(h_grp, "NumFilesPerSnapshot", 1);
io_write_attribute_i(h_grp, "ThisFile", e->nodeID); io_write_attribute_i(h_grp, "ThisFile", 0);
io_write_attribute_s(h_grp, "SelectOutput", "Default"); io_write_attribute_s(h_grp, "SelectOutput", "Default");
io_write_attribute_i(h_grp, "Virtual", virtual_file); io_write_attribute_i(h_grp, "Virtual", virtual_file);
const int to_write[swift_type_count] = {0}; const int to_write[swift_type_count] = {0};
...@@ -138,219 +135,6 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e, ...@@ -138,219 +135,6 @@ void write_fof_hdf5_header(hid_t h_file, const struct engine* e,
/*fof=*/1); /*fof=*/1);
} }
void write_virtual_fof_hdf5_array(
const struct engine* e, hid_t grp, const char* fileName_base,
const char* partTypeGroupName, const struct io_props props,
const size_t N_total, const long long* N_counts,
const enum lossy_compression_schemes lossy_compression,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units) {
#if H5_VERSION_GE(1, 10, 0)
/* Create data space */
hid_t h_space;
if (N_total > 0)
h_space = H5Screate(H5S_SIMPLE);
else
h_space = H5Screate(H5S_NULL);
if (h_space < 0)
error("Error while creating data space for field '%s'.", props.name);
int rank = 0;
hsize_t shape[2];
hsize_t source_shape[2];
hsize_t start[2] = {0, 0};
hsize_t count[2];
if (props.dimension > 1) {
rank = 2;
shape[0] = N_total;
shape[1] = props.dimension;
source_shape[0] = 0;
source_shape[1] = props.dimension;
count[0] = 0;
count[1] = props.dimension;
} else {
rank = 1;
shape[0] = N_total;
shape[1] = 0;
source_shape[0] = 0;
source_shape[1] = 0;
count[0] = 0;
count[1] = 0;
}
/* Change shape of data space */
hid_t h_err = H5Sset_extent_simple(h_space, rank, shape, NULL);
if (h_err < 0)
error("Error while changing data space shape for field '%s'.", props.name);
/* Dataset type */
hid_t h_type = H5Tcopy(io_hdf5_type(props.type));
/* Dataset properties */
hid_t h_prop = H5Pcreate(H5P_DATASET_CREATE);
/* Create filters and set compression level if we have something to write */
char comp_buffer[32] = "None";
/* The name of the dataset to map to in the other files */
char source_dataset_name[256];
sprintf(source_dataset_name, "Groups/%s", props.name);
/* Construct a relative base name */
char fileName_relative_base[256];
int pos_last_slash = strlen(fileName_base) - 1;
for (/* */; pos_last_slash >= 0; --pos_last_slash)
if (fileName_base[pos_last_slash] == '/') break;
sprintf(fileName_relative_base, "%s", &fileName_base[pos_last_slash + 1]);
/* Create all the virtual mappings */
for (int i = 0; i < e->nr_nodes; ++i) {
/* Get the number of particles of this type written on this rank */
count[0] = N_counts[i];
/* Select the space in the virtual file */
h_err = H5Sselect_hyperslab(h_space, H5S_SELECT_SET, start, /*stride=*/NULL,
count, /*block=*/NULL);
if (h_err < 0) error("Error selecting hyper-slab in the virtual file");
/* Select the space in the (already existing) source file */
source_shape[0] = count[0];
hid_t h_source_space = H5Screate_simple(rank, source_shape, NULL);
if (h_source_space < 0) error("Error creating space in the source file");
char fileName[1024];
sprintf(fileName, "%s.%d.hdf5", fileName_relative_base, i);
/* Make the virtual link */
h_err = H5Pset_virtual(h_prop, h_space, fileName, source_dataset_name,
h_source_space);
if (h_err < 0) error("Error setting the virtual properties");
H5Sclose(h_source_space);
/* Move to the next slab (i.e. next file) */
start[0] += count[0];
}
/* Create virtual dataset */
const hid_t h_data = H5Dcreate(grp, props.name, h_type, h_space, H5P_DEFAULT,
h_prop, H5P_DEFAULT);
if (h_data < 0) error("Error while creating dataspace '%s'.", props.name);
/* Write unit conversion factors for this data set */
char buffer[FIELD_BUFFER_SIZE] = {0};
units_cgs_conversion_string(buffer, snapshot_units, props.units,
props.scale_factor_exponent);
float baseUnitsExp[5];
units_get_base_unit_exponents_array(baseUnitsExp, props.units);
io_write_attribute_f(h_data, "U_M exponent", baseUnitsExp[UNIT_MASS]);
io_write_attribute_f(h_data, "U_L exponent", baseUnitsExp[UNIT_LENGTH]);
io_write_attribute_f(h_data, "U_t exponent", baseUnitsExp[UNIT_TIME]);
io_write_attribute_f(h_data, "U_I exponent", baseUnitsExp[UNIT_CURRENT]);
io_write_attribute_f(h_data, "U_T exponent", baseUnitsExp[UNIT_TEMPERATURE]);
io_write_attribute_f(h_data, "h-scale exponent", 0.f);
io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
/* Write the actual number this conversion factor corresponds to */
const double factor =
units_cgs_conversion_factor(snapshot_units, props.units);
io_write_attribute_d(
h_data,
"Conversion factor to CGS (not including cosmological corrections)",
factor);
io_write_attribute_d(
h_data,
"Conversion factor to physical CGS (including cosmological corrections)",
factor * pow(e->cosmology->a, props.scale_factor_exponent));
#ifdef SWIFT_DEBUG_CHECKS
if (strlen(props.description) == 0)
error("Invalid (empty) description of the field '%s'", props.name);
#endif
/* Write the full description */
io_write_attribute_s(h_data, "Description", props.description);
/* Close everything */
H5Tclose(h_type);
H5Pclose(h_prop);
H5Dclose(h_data);
H5Sclose(h_space);
#endif
}
void write_fof_virtual_file(const struct fof_props* props,
const size_t num_groups_total,
const long long* N_counts, const struct engine* e) {
#if H5_VERSION_GE(1, 10, 0)
/* Create the filename */
char file_name[512];
char file_name_base[512];
char subdir_name[265];
sprintf(subdir_name, "%s_%04i", props->base_name, e->snapshot_output_count);
const char* base = basename(subdir_name);
sprintf(file_name_base, "%s/%s", subdir_name, base);
sprintf(file_name, "%s/%s.hdf5", subdir_name, base);
/* Open HDF5 file with the chosen parameters */
hid_t h_file = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
if (h_file < 0) error("Error while opening file '%s'.", file_name);
/* Start by writing the header */
write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_total, props,
/*virtual_file=*/1);
hid_t h_grp =
H5Gcreate(h_file, "/Groups", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
if (h_grp < 0) error("Error while creating groups group.\n");
struct io_props output_prop;
output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f,
(char*)props->group_mass, sizeof(double),
"FOF group masses");
write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
num_groups_total, N_counts,
compression_write_lossless, e->internal_units,
e->snapshot_units);
output_prop =
io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f,
(char*)props->group_centre_of_mass,
3 * sizeof(double), "FOF group centres of mass");
write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
num_groups_total, N_counts,
compression_write_lossless, e->internal_units,
e->snapshot_units);
output_prop = io_make_output_field_(
"GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
(char*)props->group_index, sizeof(size_t), "FOF group IDs");
write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
num_groups_total, N_counts,
compression_write_lossless, e->internal_units,
e->snapshot_units);
output_prop =
io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
(char*)props->final_group_size, sizeof(long long),
"FOF group length (number of particles)");
write_virtual_fof_hdf5_array(e, h_grp, file_name_base, "Groups", output_prop,
num_groups_total, N_counts,
compression_write_lossless, e->internal_units,
e->snapshot_units);
/* Close everything */
H5Gclose(h_grp);
H5Fclose(h_file);
#endif
}
void write_fof_hdf5_array( void write_fof_hdf5_array(
const struct engine* e, hid_t grp, const char* fileName, const struct engine* e, hid_t grp, const char* fileName,
const char* partTypeGroupName, const struct io_props props, const size_t N, const char* partTypeGroupName, const struct io_props props, const size_t N,
...@@ -476,6 +260,9 @@ void write_fof_hdf5_array( ...@@ -476,6 +260,9 @@ void write_fof_hdf5_array(
io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent); io_write_attribute_f(h_data, "a-scale exponent", props.scale_factor_exponent);
io_write_attribute_s(h_data, "Expression for physical CGS units", buffer); io_write_attribute_s(h_data, "Expression for physical CGS units", buffer);
io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer); io_write_attribute_s(h_data, "Lossy compression filter", comp_buffer);
io_write_attribute_b(h_data, "Value stored as physical", props.is_physical);
io_write_attribute_b(h_data, "Property can be converted to comoving",
props.is_convertible_to_comoving);
/* Write the actual number this conversion factor corresponds to */ /* Write the actual number this conversion factor corresponds to */
const double factor = const double factor =
...@@ -506,37 +293,24 @@ void write_fof_hdf5_array( ...@@ -506,37 +293,24 @@ void write_fof_hdf5_array(
} }
void write_fof_hdf5_catalogue(const struct fof_props* props, void write_fof_hdf5_catalogue(const struct fof_props* props,
const size_t num_groups, const struct engine* e) { const struct engine* e) {
char file_name[512]; char file_name[512];
#ifdef WITH_MPI
char subdir_name[265];
sprintf(subdir_name, "%s_%04i", props->base_name, e->snapshot_output_count);
if (e->nodeID == 0) safe_checkdir(subdir_name, /*create=*/1);
MPI_Barrier(MPI_COMM_WORLD);
const char* base = basename(subdir_name);
sprintf(file_name, "%s/%s.%d.hdf5", subdir_name, base, e->nodeID);
#else
sprintf(file_name, "%s_%04i.hdf5", props->base_name, sprintf(file_name, "%s_%04i.hdf5", props->base_name,
e->snapshot_output_count); e->snapshot_output_count);
#endif
hid_t h_file = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); /* Set the minimal API version to avoid issues with advanced features */
hid_t h_props = H5Pcreate(H5P_FILE_ACCESS);
herr_t err = H5Pset_libver_bounds(h_props, HDF5_LOWEST_FILE_FORMAT_VERSION,
HDF5_HIGHEST_FILE_FORMAT_VERSION);
if (err < 0) error("Error setting the hdf5 API version");
hid_t h_file = H5Fcreate(file_name, H5F_ACC_TRUNC, H5P_DEFAULT, h_props);
if (h_file < 0) error("Error while opening file '%s'.", file_name); if (h_file < 0) error("Error while opening file '%s'.", file_name);
/* Compute the number of groups */ /* Compute the number of groups */
long long num_groups_local = num_groups; long long num_groups_local = props->num_groups;
long long num_groups_total = num_groups; long long num_groups_total = props->num_groups;
#ifdef WITH_MPI
MPI_Allreduce(&num_groups, &num_groups_total, 1, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
/* Rank 0 collects the number of groups written by each rank */
long long* N_counts = (long long*)malloc(e->nr_nodes * sizeof(long long));
MPI_Gather(&num_groups_local, 1, MPI_LONG_LONG_INT, N_counts, 1,
MPI_LONG_LONG_INT, 0, MPI_COMM_WORLD);
#endif
/* Start by writing the header */ /* Start by writing the header */
write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_local, props, write_fof_hdf5_header(h_file, e, num_groups_total, num_groups_local, props,
...@@ -549,27 +323,32 @@ void write_fof_hdf5_catalogue(const struct fof_props* props, ...@@ -549,27 +323,32 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
struct io_props output_prop; struct io_props output_prop;
output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f, output_prop = io_make_output_field_("Masses", DOUBLE, 1, UNIT_CONV_MASS, 0.f,
(char*)props->group_mass, sizeof(double), (char*)props->group_mass, sizeof(double),
"FOF group masses"); "FOF group masses", /*physical=*/0,
/*convertible_to_comoving=*/1);
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop, write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless, num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units); e->internal_units, e->snapshot_units);
output_prop = output_prop =
io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, io_make_output_field_("Centres", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f,
(char*)props->group_centre_of_mass, (char*)props->group_centre_of_mass,
3 * sizeof(double), "FOF group centres of mass"); 3 * sizeof(double), "FOF group centres of mass",
/*physical=*/0, /*convertible_to_comoving=*/1);
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop, write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless, num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units); e->internal_units, e->snapshot_units);
output_prop = io_make_output_field_( output_prop =
"GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, io_make_output_field_("GroupIDs", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
(char*)props->group_index, sizeof(size_t), "FOF group IDs"); (char*)props->final_group_index, sizeof(long long),
"FOF group IDs", /*physical=*/1,
/*convertible_to_comoving=*/0);
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop, write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless, num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units); e->internal_units, e->snapshot_units);
output_prop = output_prop =
io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f, io_make_output_field_("Sizes", LONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
(char*)props->final_group_size, sizeof(long long), (char*)props->final_group_size, sizeof(long long),
"FOF group length (number of particles)"); "FOF group length (number of particles)",
/*physical=*/1, /*convertible_to_comoving=*/0);
write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop, write_fof_hdf5_array(e, h_grp, file_name, "Groups", output_prop,
num_groups_local, compression_write_lossless, num_groups_local, compression_write_lossless,
e->internal_units, e->snapshot_units); e->internal_units, e->snapshot_units);
...@@ -577,20 +356,7 @@ void write_fof_hdf5_catalogue(const struct fof_props* props, ...@@ -577,20 +356,7 @@ void write_fof_hdf5_catalogue(const struct fof_props* props,
/* Close everything */ /* Close everything */
H5Gclose(h_grp); H5Gclose(h_grp);
H5Fclose(h_file); H5Fclose(h_file);
H5Pclose(h_props);
#ifdef WITH_MPI
#if H5_VERSION_GE(1, 10, 0)
/* Write the virtual meta-file */
if (e->nodeID == 0)
write_fof_virtual_file(props, num_groups_total, N_counts, e);
#endif
/* Free the counts-per-rank array */
free(N_counts);
MPI_Barrier(MPI_COMM_WORLD);
#endif
} }
#endif /* HAVE_HDF5 */ #endif /* HAVE_HDF5 */
...@@ -25,7 +25,7 @@ ...@@ -25,7 +25,7 @@
#ifdef WITH_FOF #ifdef WITH_FOF
void write_fof_hdf5_catalogue(const struct fof_props *props, void write_fof_hdf5_catalogue(const struct fof_props *props,
const size_t num_groups, const struct engine *e); const struct engine *e);
#endif /* WITH_FOF */ #endif /* WITH_FOF */
......
...@@ -34,6 +34,8 @@ ...@@ -34,6 +34,8 @@
#include "./forcing/roberts_flow/forcing.h" #include "./forcing/roberts_flow/forcing.h"
#elif defined(FORCING_ROBERTS_FLOW_ACCELERATION) #elif defined(FORCING_ROBERTS_FLOW_ACCELERATION)
#include "./forcing/roberts_flow_acceleration/forcing.h" #include "./forcing/roberts_flow_acceleration/forcing.h"
#elif defined(FORCING_ABC_FLOW)
#include "./forcing/ABC_flow/forcing.h"
#else #else
#error "Invalid choice of forcing terms" #error "Invalid choice of forcing terms"
#endif #endif
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2023 Matthieu Schaller (schaller@strw.leidenuniv.nl)
* Copyright (c) 2024 Nikyta Shchutksyi (shchutskyi@lorentz.leidenuniv.nl)
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_FORCING_ABC_FLOW_H
#define SWIFT_FORCING_ABC_FLOW_H
/* Config parameters. */
#include <config.h>
/* Standard includes. */
#include <float.h>
/* Local includes. */
#include "error.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "space.h"
#include "units.h"
/**
* @brief Forcing Term Properties
*/
struct forcing_terms {
/*! Reference velocity (internal units) */
float u0;
/*! Velocity scaling along the z direction */
float Vz_factor;
/*! Wavenumber of the flow */
float kv;
/*! ABC flow coefficients */
float A, B, C;
};
/**
* @brief Computes the forcing terms.
*
* Based on David Galloway (2012) ABC flows then and now, Geophysical &
* Astrophysical Fluid Dynamics, 106:4-5, 450-467 This version differs from the
* paper by imposing normalized velocity
*
* @param time The current time.
* @param terms The properties of the forcing terms.
* @param s The #space we act on.
* @param phys_const The physical constants in internal units.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void forcing_terms_apply(
const double time, const struct forcing_terms* terms, const struct space* s,
const struct phys_const* phys_const, struct part* p, struct xpart* xp) {
const double L = s->dim[0];
const float u0 = terms->u0;
const float A = terms->A;
const float B = terms->B;
const float C = terms->C;
const float Norm = 1 / sqrt(A * A + B * B + C * C);
const float Vz_factor = terms->Vz_factor;
const double k0 = (2. * M_PI / L) * terms->kv;
double v_ABC[3];
/* Eq. 2 of David Galloway (2012) ABC flows then and now, Geophysical &
* Astrophysical Fluid Dynamics, 106:4-5, 450-467 */
// Velocity normalized such that <v>rms = u0
v_ABC[0] = u0 * Norm * (A * sin(k0 * p->x[2]) + C * cos(k0 * p->x[1]));
v_ABC[1] = u0 * Norm * (B * sin(k0 * p->x[0]) + A * cos(k0 * p->x[2]));
v_ABC[2] = u0 * Norm * (C * sin(k0 * p->x[1]) + B * cos(k0 * p->x[0]));
/* Force the velocity and possibly scale the z-direction */
xp->v_full[0] = v_ABC[0];
xp->v_full[1] = v_ABC[1];
xp->v_full[2] = v_ABC[2] * Vz_factor;
}
/**
* @brief Computes the time-step condition due to the forcing terms.
*
* Nothing to do here. --> Return FLT_MAX.
*
* @param time The current time.
* @param terms The properties of the forcing terms.
* @param phys_const The physical constants in internal units.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static float forcing_terms_timestep(
double time, const struct forcing_terms* terms,
const struct phys_const* phys_const, const struct part* p,
const struct xpart* xp) {
return FLT_MAX;
}
/**
* @brief Prints the properties of the forcing terms to stdout.
*
* @param terms The #forcing_terms properties of the run.
*/
static INLINE void forcing_terms_print(const struct forcing_terms* terms) {
message("Forcing terms is 'ABC flow'. U0: %.5f / Vz factor: %.5f.", terms->u0,
terms->Vz_factor);
message("Run using ABC parameters: A = %.5f, B = %.5f, C = %.5f", terms->A,
terms->B, terms->C);
}
/**
* @brief Initialises the forcing term properties
*
* @param parameter_file The parsed parameter file
* @param phys_const Physical constants in internal units
* @param us The current internal system of units
* @param s The #space object.
* @param terms The forcing term properties to initialize
*/
static INLINE void forcing_terms_init(struct swift_params* parameter_file,
const struct phys_const* phys_const,
const struct unit_system* us,
const struct space* s,
struct forcing_terms* terms) {
terms->u0 = parser_get_param_double(parameter_file, "ABC_Flow_Forcing:u0");
terms->Vz_factor = parser_get_opt_param_float(
parameter_file, "ABC_Flow_Forcing:Vz_factor", 1.f);
terms->kv = parser_get_param_double(parameter_file, "ABC_Flow_Forcing:kv");
terms->A = parser_get_param_double(parameter_file, "ABC_Flow_Forcing:A");
terms->B = parser_get_param_double(parameter_file, "ABC_Flow_Forcing:B");
terms->C = parser_get_param_double(parameter_file, "ABC_Flow_Forcing:C");
}
#endif /* SWIFT_FORCING_ABC_FLOW_H */