diff --git a/examples/main.c b/examples/main.c index f5622b0615871d72ffadb8dedb50b0de4271fa31..6c72e6c66c9cbf08364ca5df2900926bea9bf8cd 100644 --- a/examples/main.c +++ b/examples/main.c @@ -735,11 +735,10 @@ int main(int argc, char *argv[]) { error("ERROR: Running with feedback but compiled without it."); #endif feedback_props_init(&feedback_properties, &prog_const, &us, params, - &hydro_properties, &cosmo); - } - else + &hydro_properties, &cosmo); + } else bzero(&feedback_properties, sizeof(struct feedback_props)); - + /* Initialise the gravity properties */ if (with_self_gravity) gravity_props_init(&gravity_properties, params, &cosmo, with_cosmology, @@ -970,7 +969,7 @@ int main(int argc, char *argv[]) { engine_policies, talking, &reparttype, &us, &prog_const, &cosmo, &hydro_properties, &entropy_floor, &gravity_properties, &stars_properties, &feedback_properties, &mesh, &potential, - &cooling_func, &starform, &chemistry); + &cooling_func, &starform, &chemistry); engine_config(0, &e, params, nr_nodes, myrank, nr_threads, with_aff, talking, restart_file); diff --git a/src/engine.c b/src/engine.c index a2ad8816e88f94b3b55f6d9d01e691f04b5cfc02..3be51280d9785990cea8121409087f0f935e2a4f 100644 --- a/src/engine.c +++ b/src/engine.c @@ -4404,8 +4404,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, struct cosmology *cosmo, struct hydro_props *hydro, const struct entropy_floor_properties *entropy_floor, struct gravity_props *gravity, const struct stars_props *stars, - const struct feedback_props *feedback, - struct pm_mesh *mesh, + const struct feedback_props *feedback, struct pm_mesh *mesh, const struct external_potential *potential, struct cooling_function_data *cooling_func, const struct star_formation *starform, diff --git a/src/engine.h b/src/engine.h index 88f2ae4f9d510c9701ae16b4a769982c58ceb522..fbbf82d339241dbf87a28b1ea847fc8e843bf8d5 100644 --- a/src/engine.h +++ b/src/engine.h @@ -390,7 +390,7 @@ struct engine { /* Properties of the sellar feedback model */ const struct feedback_props *feedback_props; - + /* Properties of the chemistry model */ const struct chemistry_global_data *chemistry; @@ -447,8 +447,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params, struct cosmology *cosmo, struct hydro_props *hydro, const struct entropy_floor_properties *entropy_floor, struct gravity_props *gravity, const struct stars_props *stars, - const struct feedback_props *feedback, - struct pm_mesh *mesh, + const struct feedback_props *feedback, struct pm_mesh *mesh, const struct external_potential *potential, struct cooling_function_data *cooling_func, const struct star_formation *starform, diff --git a/src/feedback.h b/src/feedback.h index e4987f916822fd90affd27372a51b6f11cf73287..d988e0280ee370947277ba909244161757743586 100644 --- a/src/feedback.h +++ b/src/feedback.h @@ -24,8 +24,8 @@ /* Select the correct feedback model */ #if defined(FEEDBACK_NONE) -#include "./feedback/Default/feedback.h" #include "./deedback/Default/feedback_iact.h" +#include "./feedback/Default/feedback.h" #elif defined(FEEDBACK_EAGLE) #include "./feedback/EAGLE/feedback.h" #include "./feedback/EAGLE/feedback_iact.h" diff --git a/src/feedback/EAGLE/feedback.c b/src/feedback/EAGLE/feedback.c index 3cf1313c52d1a68e787727e2126935b26f9a5a04..08a75e55278c7c1c321ced809b99adec7decb9a7 100644 --- a/src/feedback/EAGLE/feedback.c +++ b/src/feedback/EAGLE/feedback.c @@ -17,7 +17,6 @@ * ******************************************************************************/ - #include "feedback.h" #include "hydro_properties.h" @@ -34,14 +33,14 @@ * @param dt length of step */ float compute_SNe(const struct spart* sp, - const struct feedback_props* feedback_props, - const float age, const float dt) { + const struct feedback_props* feedback_props, const float age, + const float dt) { const float SNII_wind_delay = feedback_props->SNII_wind_delay; - + if (age <= SNII_wind_delay && (age + dt) > SNII_wind_delay) { - return feedback_props->num_SNII_per_msun * sp->mass_init / + return feedback_props->num_SNII_per_msun * sp->mass_init / feedback_props->const_solar_mass; } else { return 0.f; @@ -65,13 +64,14 @@ inline static void determine_bin_yield_AGB( const struct feedback_props* feedback_props) { // MATTHIEU - + /* if (log10_metallicity > log10_min_metallicity) { */ /* /\* Find metallicity bin which contains the star's metallicity *\/ */ /* int j; */ /* for (j = 0; j < star_properties->feedback.AGB_n_z - 1 && */ /* log10_metallicity > */ - /* star_properties->feedback.yield_AGB.metallicity[j + 1]; */ + /* star_properties->feedback.yield_AGB.metallicity[j + 1]; + */ /* j++) */ /* ; */ /* *iz_low = j; */ @@ -89,8 +89,10 @@ inline static void determine_bin_yield_AGB( /* *dz = 0; */ /* /\* Normalize offset *\/ */ - /* float deltaz = star_properties->feedback.yield_AGB.metallicity[*iz_high] - */ - /* star_properties->feedback.yield_AGB.metallicity[*iz_low]; */ + /* float deltaz = star_properties->feedback.yield_AGB.metallicity[*iz_high] + * - */ + /* star_properties->feedback.yield_AGB.metallicity[*iz_low]; + */ /* if (deltaz > 0) */ /* *dz /= deltaz; */ @@ -120,13 +122,14 @@ inline static void determine_bin_yield_SNII( const struct feedback_props* feedback_props) { // MATTHIEU - + /* if (log10_metallicity > log10_min_metallicity) { */ /* /\* Find metallicity bin which contains the star's metallicity *\/ */ /* int j; */ /* for (j = 0; j < star_properties->feedback.SNII_n_z - 1 && */ /* log10_metallicity > */ - /* star_properties->feedback.yield_SNII.metallicity[j + 1]; */ + /* star_properties->feedback.yield_SNII.metallicity[j + 1]; + */ /* j++) */ /* ; */ /* *iz_low = j; */ @@ -144,8 +147,10 @@ inline static void determine_bin_yield_SNII( /* *dz = 0; */ /* /\* Normalize offset *\/ */ - /* float deltaz = star_properties->feedback.yield_SNII.metallicity[*iz_high] - */ - /* star_properties->feedback.yield_SNII.metallicity[*iz_low]; */ + /* float deltaz = star_properties->feedback.yield_SNII.metallicity[*iz_high] + * - */ + /* star_properties->feedback.yield_SNII.metallicity[*iz_low]; + */ /* if (deltaz > 0) */ /* *dz = *dz / deltaz; */ @@ -158,7 +163,6 @@ inline static void determine_bin_yield_SNII( /* } */ } - /** * @brief compute enrichment and feedback due to SNIa. To do this compute the * number of SNIa that occur during the timestep, multiply by constants read @@ -183,9 +187,9 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass, * use updated values for the star's age and timestep in this function */ if (log10_max_mass > feedback_props->log10_SNIa_max_mass_msun) { log10_max_mass = feedback_props->log10_SNIa_max_mass_msun; - float lifetime_Gyr = - lifetime_in_Gyr(exp(M_LN10 * feedback_props->log10_SNIa_max_mass_msun), - sp->chemistry_data.metal_mass_fraction_total, feedback_props); + float lifetime_Gyr = lifetime_in_Gyr( + exp(M_LN10 * feedback_props->log10_SNIa_max_mass_msun), + sp->chemistry_data.metal_mass_fraction_total, feedback_props); dt_Gyr = star_age_Gyr + dt_Gyr - lifetime_Gyr; star_age_Gyr = lifetime_Gyr; } @@ -196,7 +200,7 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass, feedback_props->SNIa_efficiency * (exp(-star_age_Gyr / feedback_props->SNIa_timescale_Gyr) - exp(-(star_age_Gyr + dt_Gyr) / feedback_props->SNIa_timescale_Gyr)) * - sp->mass_init; + sp->mass_init; sp->feedback_data.to_distribute.num_SNIa = num_SNIa_per_msun / feedback_props->const_solar_mass; @@ -237,11 +241,11 @@ inline static void evolve_SNIa(float log10_min_mass, float log10_max_mass, * @param stars star properties data structure * @param sp spart we are computing feedback from */ -inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, - float* stellar_yields, - const struct feedback_props* restrict feedback_props, - struct spart* restrict sp) { - +inline static void evolve_SNII( + float log10_min_mass, float log10_max_mass, float* stellar_yields, + const struct feedback_props* restrict feedback_props, + struct spart* restrict sp) { + int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index; /* If mass at beginning of step is less than tabulated lower bound for IMF, @@ -263,11 +267,12 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, &high_imf_mass_bin_index, feedback_props); /* Integrate IMF to determine number of SNII */ - sp->feedback_data.to_distribute.num_SNII = - integrate_imf(log10_min_mass, log10_max_mass, 0, stellar_yields, feedback_props); + sp->feedback_data.to_distribute.num_SNII = integrate_imf( + log10_min_mass, log10_max_mass, 0, stellar_yields, feedback_props); /* determine which metallicity bin and offset this star belongs to */ - int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, high_index_2d; + int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, + high_index_2d; float dz = 0.; determine_bin_yield_SNII(&iz_low, &iz_high, &dz, log10(sp->chemistry_data.metal_mass_fraction_total), @@ -284,12 +289,12 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, high_index_3d = row_major_index_3d( iz_high, elem, mass_bin_index, feedback_props->SNII_n_z, chemistry_element_count, feedback_props->n_imf_mass_bins); - low_index_2d = row_major_index_2d( - iz_low, mass_bin_index, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * (feedback_props->yield_SNII.yield_IMF_resampled[low_index_3d] + @@ -302,19 +307,19 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, .ejecta_IMF_resampled[high_index_2d]); } - metal_mass_released[elem] = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + metal_mass_released[elem] = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); } /* Compute mass produced */ for (mass_bin_index = low_imf_mass_bin_index; mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) { - low_index_2d = row_major_index_2d(iz_low, mass_bin_index, - feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * (feedback_props->yield_SNII .total_metals_IMF_resampled[low_index_2d] + @@ -328,8 +333,8 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, .ejecta_IMF_resampled[high_index_2d]); } - metal_mass_released_total = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + metal_mass_released_total = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); /* yield normalization */ float mass_ejected, mass_released; @@ -343,20 +348,20 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, /* compute the total metal mass ejected from the star*/ for (mass_bin_index = low_imf_mass_bin_index; mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) { - low_index_2d = row_major_index_2d(iz_low, mass_bin_index, - feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * feedback_props->yield_SNII.ejecta_IMF_resampled[low_index_2d] + dz * feedback_props->yield_SNII.ejecta_IMF_resampled[high_index_2d]; } - mass_ejected = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + mass_ejected = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); /* compute the total mass released */ mass_released = metal_mass_released_total + @@ -370,10 +375,12 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, const float norm_factor = sp->mass_init * mass_ejected / mass_released; for (int i = 0; i < chemistry_element_count; i++) { - sp->feedback_data.to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor; + sp->feedback_data.to_distribute.metal_mass[i] += + metal_mass_released[i] * norm_factor; } for (int i = 0; i < chemistry_element_count; i++) { - sp->feedback_data.to_distribute.mass_from_SNII += sp->feedback_data.to_distribute.metal_mass[i]; + sp->feedback_data.to_distribute.mass_from_SNII += + sp->feedback_data.to_distribute.metal_mass[i]; } sp->feedback_data.to_distribute.total_metal_mass += metal_mass_released_total * norm_factor; @@ -396,11 +403,11 @@ inline static void evolve_SNII(float log10_min_mass, float log10_max_mass, * @param feedback_props star properties data structure * @param sp spart we are computing feedback from */ -inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, - float* stellar_yields, - const struct feedback_props* restrict feedback_props, - struct spart* restrict sp) { - +inline static void evolve_AGB( + float log10_min_mass, float log10_max_mass, float* stellar_yields, + const struct feedback_props* restrict feedback_props, + struct spart* restrict sp) { + int low_imf_mass_bin_index, high_imf_mass_bin_index, mass_bin_index; /* If mass at end of step is greater than tabulated lower bound for IMF, limit @@ -417,7 +424,8 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, &high_imf_mass_bin_index, feedback_props); /* determine which metallicity bin and offset this star belongs to */ - int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, high_index_2d; + int iz_low = 0, iz_high = 0, low_index_3d, high_index_3d, low_index_2d, + high_index_2d; float dz = 0.f; determine_bin_yield_AGB(&iz_low, &iz_high, &dz, log10(sp->chemistry_data.metal_mass_fraction_total), @@ -434,12 +442,12 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, high_index_3d = row_major_index_3d( iz_high, elem, mass_bin_index, feedback_props->AGB_n_z, chemistry_element_count, feedback_props->n_imf_mass_bins); - low_index_2d = row_major_index_2d( - iz_low, mass_bin_index, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * (feedback_props->yield_AGB.yield_IMF_resampled[low_index_3d] + @@ -452,19 +460,19 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, .ejecta_IMF_resampled[high_index_2d]); } - metal_mass_released[elem] = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + metal_mass_released[elem] = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); } /* Compute mass produced */ for (mass_bin_index = low_imf_mass_bin_index; mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) { - low_index_2d = row_major_index_2d(iz_low, mass_bin_index, - feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * (feedback_props->yield_AGB @@ -478,8 +486,8 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, feedback_props->yield_AGB.ejecta_IMF_resampled[high_index_2d]); } - metal_mass_released_total = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + metal_mass_released_total = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); /* yield normalization */ float mass_ejected, mass_released; @@ -493,20 +501,20 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, /* compute the total metal mass ejected from the star */ for (mass_bin_index = low_imf_mass_bin_index; mass_bin_index < high_imf_mass_bin_index + 1; mass_bin_index++) { - low_index_2d = row_major_index_2d(iz_low, mass_bin_index, - feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); - high_index_2d = row_major_index_2d( - iz_high, mass_bin_index, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); + low_index_2d = + row_major_index_2d(iz_low, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); + high_index_2d = + row_major_index_2d(iz_high, mass_bin_index, feedback_props->AGB_n_z, + feedback_props->n_imf_mass_bins); stellar_yields[mass_bin_index] = (1 - dz) * feedback_props->yield_AGB.ejecta_IMF_resampled[low_index_2d] + dz * feedback_props->yield_AGB.ejecta_IMF_resampled[high_index_2d]; } - mass_ejected = - integrate_imf(log10_min_mass, log10_max_mass, 2, stellar_yields, feedback_props); + mass_ejected = integrate_imf(log10_min_mass, log10_max_mass, 2, + stellar_yields, feedback_props); /* compute the total mass released */ mass_released = metal_mass_released_total + @@ -520,8 +528,10 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, const float norm_factor = sp->mass_init * mass_ejected / mass_released; for (int i = 0; i < chemistry_element_count; i++) { - sp->feedback_data.to_distribute.metal_mass[i] += metal_mass_released[i] * norm_factor; - sp->feedback_data.to_distribute.mass_from_AGB += metal_mass_released[i] * norm_factor; + sp->feedback_data.to_distribute.metal_mass[i] += + metal_mass_released[i] * norm_factor; + sp->feedback_data.to_distribute.mass_from_AGB += + metal_mass_released[i] * norm_factor; } sp->feedback_data.to_distribute.total_metal_mass += metal_mass_released_total * norm_factor; @@ -543,14 +553,13 @@ inline static void evolve_AGB(float log10_min_mass, float log10_max_mass, * @param dt length of current timestep */ void compute_stellar_evolution(const struct feedback_props* feedback_props, - const struct cosmology *cosmo, struct spart* sp, - const struct unit_system* us, const float age, - const float dt) { + const struct cosmology* cosmo, struct spart* sp, + const struct unit_system* us, const float age, + const float dt) { /* Allocate temporary array for calculating imf weights */ float* stellar_yields; - stellar_yields = - malloc(feedback_props->n_imf_mass_bins * sizeof(float)); + stellar_yields = malloc(feedback_props->n_imf_mass_bins * sizeof(float)); /* Convert dt and stellar age from internal units to Gyr. */ const double Gyr_in_cgs = 1e9 * 365. * 24. * 3600.; @@ -561,11 +570,13 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, /* Get the metallicity */ const float Z = sp->chemistry_data.metal_mass_fraction_total; - + /* calculate mass of stars that has died from the star's birth up to the * beginning and end of timestep */ - const float max_dying_mass_Msun = dying_mass_msun(star_age_Gyr, Z, feedback_props); - const float min_dying_mass_Msun = dying_mass_msun(star_age_Gyr + dt_Gyr, Z, feedback_props); + const float max_dying_mass_Msun = + dying_mass_msun(star_age_Gyr, Z, feedback_props); + const float min_dying_mass_Msun = + dying_mass_msun(star_age_Gyr + dt_Gyr, Z, feedback_props); #ifdef SWIFT_DEBUG_CHECK /* Sanity check. Worth investigating if necessary as functions for evaluating @@ -573,7 +584,7 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, if (min_dying_mass_Msun > max_dying_mass_Msun) error("min dying mass is greater than max dying mass"); #endif - + /* Integration interval is zero - this can happen if minimum and maximum * dying masses are above imf_max_mass_Msun. Return without doing any * feedback. */ @@ -582,8 +593,8 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, /* Life is better in log */ const float log10_max_dying_mass_Msun = log10f(max_dying_mass_Msun); const float log10_min_dying_mass_Msun = log10f(min_dying_mass_Msun); - - /* Compute elements, energy and momentum to distribute from the + + /* Compute elements, energy and momentum to distribute from the * three channels SNIa, SNII, AGB */ evolve_SNIa(log10_min_dying_mass_Msun, log10_max_dying_mass_Msun, feedback_props, sp, star_age_Gyr, dt_Gyr); @@ -593,29 +604,30 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, stellar_yields, feedback_props, sp); /* Compute the total mass to distribute (H + He metals) */ - sp->feedback_data.to_distribute.mass = sp->feedback_data.to_distribute.total_metal_mass + - sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] + - sp->feedback_data.to_distribute.metal_mass[chemistry_element_He]; + sp->feedback_data.to_distribute.mass = + sp->feedback_data.to_distribute.total_metal_mass + + sp->feedback_data.to_distribute.metal_mass[chemistry_element_H] + + sp->feedback_data.to_distribute.metal_mass[chemistry_element_He]; /* Compute the number of type II SNe that went off */ - sp->feedback_data.to_distribute.num_SNe = compute_SNe(sp, feedback_props, age, dt); + sp->feedback_data.to_distribute.num_SNe = + compute_SNe(sp, feedback_props, age, dt); /* Compute energy change due to thermal and kinetic energy of ejecta */ sp->feedback_data.to_distribute.d_energy = sp->feedback_data.to_distribute.mass * - (feedback_props->ejecta_specific_thermal_energy + - 0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) * - cosmo->a2_inv); + (feedback_props->ejecta_specific_thermal_energy + + 0.5 * (sp->v[0] * sp->v[0] + sp->v[1] * sp->v[1] + sp->v[2] * sp->v[2]) * + cosmo->a2_inv); /* Compute probability of heating neighbouring particles */ if (dt > 0 && sp->feedback_data.ngb_mass > 0) sp->feedback_data.to_distribute.heating_probability = feedback_props->total_energy_SNe * sp->feedback_data.to_distribute.num_SNe / - (feedback_props->temp_to_u_factor * - feedback_props->SNe_deltaT_desired * sp->feedback_data.ngb_mass); + (feedback_props->temp_to_u_factor * feedback_props->SNe_deltaT_desired * + sp->feedback_data.ngb_mass); - /* Clean up */ free(stellar_yields); } @@ -631,12 +643,12 @@ void compute_stellar_evolution(const struct feedback_props* feedback_props, * @param params The parsed parameters. * @param p The already read-in properties of the hydro scheme. */ -void feedback_props_init(struct feedback_props *fp, - const struct phys_const *phys_const, - const struct unit_system *us, - struct swift_params *params, - const struct hydro_props *hydro_props, - const struct cosmology *cosmo) { +void feedback_props_init(struct feedback_props* fp, + const struct phys_const* phys_const, + const struct unit_system* us, + struct swift_params* params, + const struct hydro_props* hydro_props, + const struct cosmology* cosmo) { /* Read SNIa timscale */ fp->SNIa_timescale_Gyr = @@ -675,7 +687,8 @@ void feedback_props_init(struct feedback_props *fp, /* Calculate temperature to internal energy conversion factor */ fp->temp_to_u_factor = phys_const->const_boltzmann_k / - (hydro_props->mu_ionised * (hydro_gamma_minus_one)*phys_const->const_proton_mass); + (hydro_props->mu_ionised * + (hydro_gamma_minus_one)*phys_const->const_proton_mass); /* Read birth time to set all stars in ICs to (defaults to -1 to indicate star * present in ICs) */ @@ -685,7 +698,6 @@ void feedback_props_init(struct feedback_props *fp, /* Copy over solar mass */ fp->const_solar_mass = phys_const->const_solar_mass; - /* Set number of elements found in yield tables */ fp->SNIa_n_elements = 42; fp->SNII_n_mass = 11; @@ -706,8 +718,7 @@ void feedback_props_init(struct feedback_props *fp, /* Yield table filepath */ parser_get_param_string(params, "EAGLEFeedback:filename", fp->yield_table_path); - parser_get_param_string(params, "EAGLEFeedback:imf_model", - fp->IMF_Model); + parser_get_param_string(params, "EAGLEFeedback:imf_model", fp->IMF_Model); /* Initialise IMF */ init_imf(fp); @@ -731,8 +742,7 @@ void feedback_props_init(struct feedback_props *fp, /* Set yield_mass_bins array */ const float imf_log10_mass_bin_size = - (fp->log10_imf_max_mass_msun - - fp->log10_imf_min_mass_msun) / + (fp->log10_imf_max_mass_msun - fp->log10_imf_min_mass_msun) / (fp->n_imf_mass_bins - 1); for (int i = 0; i < fp->n_imf_mass_bins; i++) fp->yield_mass_bins[i] = @@ -748,12 +758,9 @@ void feedback_props_init(struct feedback_props *fp, /* Calculate number of type II SN per solar mass * Note: since we are integrating the IMF without weighting it by the yields * pass NULL pointer for stellar_yields array */ - fp->num_SNII_per_msun = - integrate_imf(fp->log10_SNII_min_mass_msun, - fp->log10_SNII_max_mass_msun, 0, - /*(stellar_yields=)*/ NULL, fp); + fp->num_SNII_per_msun = integrate_imf(fp->log10_SNII_min_mass_msun, + fp->log10_SNII_max_mass_msun, 0, + /*(stellar_yields=)*/ NULL, fp); message("initialized stellar feedback"); } - - diff --git a/src/feedback/EAGLE/feedback.h b/src/feedback/EAGLE/feedback.h index ea79dc96ed3581ed89267fd9754d22e0d47c7dc7..d031cfaca037555d1d794366678e9c509426ef03 100644 --- a/src/feedback/EAGLE/feedback.h +++ b/src/feedback/EAGLE/feedback.h @@ -21,16 +21,15 @@ #include "cosmology.h" #include "hydro_properties.h" -#include "units.h" #include "part.h" +#include "units.h" #include "feedback_properties.h" - void compute_stellar_evolution(const struct feedback_props* feedback_props, - const struct cosmology *cosmo, - struct spart* sp, const struct unit_system* us, - const float age, const float dt); + const struct cosmology* cosmo, struct spart* sp, + const struct unit_system* us, const float age, + const float dt); /** * @brief Prepares a s-particle for its feedback interactions @@ -114,5 +113,4 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart( sp->mass -= sp->feedback_data.to_distribute.mass; } - #endif /* SWIFT_FEEDBACK_EAGLE_H */ diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h index 3533e13b462e8eae034ca48dad2758b3ee27549a..4b11ec3cf5f9268030e451742268a9b2465d8dca 100644 --- a/src/feedback/EAGLE/feedback_iact.h +++ b/src/feedback/EAGLE/feedback_iact.h @@ -14,11 +14,13 @@ * @param ti_current Current integer time value */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_feedback_density( - float r2, const float *dx, float hi, float hj, struct spart *restrict si, - const struct part *restrict pj, const struct cosmology *restrict cosmo, - const struct feedback_props * feedback_props, - const struct xpart *restrict xp, integertime_t ti_current) { +runner_iact_nonsym_feedback_density(float r2, const float *dx, float hi, + float hj, struct spart *restrict si, + const struct part *restrict pj, + const struct cosmology *restrict cosmo, + const struct feedback_props *feedback_props, + const struct xpart *restrict xp, + integertime_t ti_current) { /* Get the gas mass. */ const float mj = hydro_get_mass(pj); @@ -41,7 +43,8 @@ runner_iact_nonsym_feedback_density( * gas particles */ const float rho = hydro_get_comoving_density(pj); - if (rho != 0) si->feedback_data.density_weighted_frac_normalisation_inv += wi / rho; + if (rho != 0) + si->feedback_data.density_weighted_frac_normalisation_inv += wi / rho; /* Compute contribution to the density */ si->rho_gas += mj * wi; @@ -93,7 +96,8 @@ runner_iact_nonsym_feedback_apply( /* Update particle mass */ const float current_mass = hydro_get_mass(pj); - const float new_mass = current_mass + si->feedback_data.to_distribute.mass * density_weighted_frac; + const float new_mass = current_mass + si->feedback_data.to_distribute.mass * + density_weighted_frac; hydro_set_mass(pj, new_mass); @@ -111,8 +115,8 @@ runner_iact_nonsym_feedback_apply( const float current_metal_mass = pj->chemistry_data.metal_mass_fraction[elem] * current_mass; const float new_metal_mass = - current_metal_mass + - si->feedback_data.to_distribute.metal_mass[elem] * density_weighted_frac; + current_metal_mass + si->feedback_data.to_distribute.metal_mass[elem] * + density_weighted_frac; pj->chemistry_data.metal_mass_fraction[elem] = new_metal_mass / new_mass; } @@ -138,7 +142,8 @@ runner_iact_nonsym_feedback_apply( pj->chemistry_data.metal_mass_fraction_from_SNIa * current_mass; const float new_metal_mass_fraction_from_SNIa = current_metal_mass_fraction_from_SNIa + - si->feedback_data.to_distribute.metal_mass_from_SNIa * density_weighted_frac; + si->feedback_data.to_distribute.metal_mass_from_SNIa * + density_weighted_frac; pj->chemistry_data.metal_mass_fraction_from_SNIa = new_metal_mass_fraction_from_SNIa / new_mass; @@ -155,7 +160,8 @@ runner_iact_nonsym_feedback_apply( pj->chemistry_data.metal_mass_fraction_from_SNII * current_mass; const float new_metal_mass_fraction_from_SNII = current_metal_mass_fraction_from_SNII + - si->feedback_data.to_distribute.metal_mass_from_SNII * density_weighted_frac; + si->feedback_data.to_distribute.metal_mass_from_SNII * + density_weighted_frac; pj->chemistry_data.metal_mass_fraction_from_SNII = new_metal_mass_fraction_from_SNII / new_mass; @@ -172,13 +178,15 @@ runner_iact_nonsym_feedback_apply( pj->chemistry_data.metal_mass_fraction_from_AGB * current_mass; const float new_metal_mass_fraction_from_AGB = current_metal_mass_fraction_from_AGB + - si->feedback_data.to_distribute.metal_mass_from_AGB * density_weighted_frac; + si->feedback_data.to_distribute.metal_mass_from_AGB * + density_weighted_frac; pj->chemistry_data.metal_mass_fraction_from_AGB = new_metal_mass_fraction_from_AGB / new_mass; /* Update momentum */ for (int i = 0; i < 3; i++) { - pj->v[i] += si->feedback_data.to_distribute.mass * density_weighted_frac * (si->v[i] - pj->v[i]); + pj->v[i] += si->feedback_data.to_distribute.mass * density_weighted_frac * + (si->v[i] - pj->v[i]); } /* Energy feedback */ @@ -189,15 +197,17 @@ runner_iact_nonsym_feedback_apply( if (feedback_props->continuous_heating) { // We're doing ONLY continuous heating d_energy += si->feedback_data.to_distribute.num_SNIa * - feedback_props->total_energy_SNe * - density_weighted_frac * si->mass_init; + feedback_props->total_energy_SNe * density_weighted_frac * + si->mass_init; } else { // We're doing stochastic heating heating_probability = si->feedback_data.to_distribute.heating_probability; du = feedback_props->SNe_deltaT_desired * feedback_props->temp_to_u_factor; - + if (heating_probability >= 1) { - du = feedback_props->total_energy_SNe * si->feedback_data.to_distribute.num_SNIa / si->feedback_data.ngb_mass; + du = feedback_props->total_energy_SNe * + si->feedback_data.to_distribute.num_SNIa / + si->feedback_data.ngb_mass; heating_probability = 1; } diff --git a/src/feedback/EAGLE/feedback_properties.h b/src/feedback/EAGLE/feedback_properties.h index b52bfddc6f3bb9ab0c50186a7a4ab903e5656e6b..66f8d02adcbe39bcb3531677bb6ba19331de33c7 100644 --- a/src/feedback/EAGLE/feedback_properties.h +++ b/src/feedback/EAGLE/feedback_properties.h @@ -150,11 +150,9 @@ struct feedback_props { float spart_first_init_birth_time; }; - - void feedback_props_init(struct feedback_props *fp, - const struct phys_const *phys_const, - const struct unit_system *us, - struct swift_params *params, - const struct hydro_props *hydro_props, - const struct cosmology *cosmo); + const struct phys_const *phys_const, + const struct unit_system *us, + struct swift_params *params, + const struct hydro_props *hydro_props, + const struct cosmology *cosmo); diff --git a/src/feedback/EAGLE/feedback_struct.h b/src/feedback/EAGLE/feedback_struct.h index 093bbfc424f507ab4cb9d670954f7cfa37a816a4..f198d851cee7e70c6c375d10f1129aff22fa8043 100644 --- a/src/feedback/EAGLE/feedback_struct.h +++ b/src/feedback/EAGLE/feedback_struct.h @@ -19,20 +19,19 @@ #ifndef SWIFT_FEEDBACK_STRUCT_EAGLE_H #define SWIFT_FEEDBACK_STRUCT_EAGLE_H - struct feedback_part_data { struct { - + /* Mass of ejecta */ float mass; - + /* Total metal mass released */ float total_metal_mass; - + /* Total mass released by element */ float metal_mass[chemistry_element_count]; - + /*! Total mass released due to SNIa */ float mass_from_SNIa; @@ -79,6 +78,4 @@ struct feedback_part_data { float ngb_mass; }; - #endif /* SWIFT_FEEDBACK_STRUCT_EAGLE_H */ - diff --git a/src/feedback/EAGLE/imf.h b/src/feedback/EAGLE/imf.h index 355a97c3de3988def133c6eff215967655aa673a..734781dd4a62c7521f1bd777fd580587ff26df98 100644 --- a/src/feedback/EAGLE/imf.h +++ b/src/feedback/EAGLE/imf.h @@ -96,8 +96,7 @@ inline static float integrate_imf(const float log10_min_mass, const int N_bins = feedback_props->n_imf_mass_bins; const float *imf = feedback_props->imf; const float *imf_mass_bin = feedback_props->imf_mass_bin; - const float *imf_mass_bin_log10 = - feedback_props->imf_mass_bin_log10; + const float *imf_mass_bin_log10 = feedback_props->imf_mass_bin_log10; /* IMF mass bin spacing in log10 space. Assumes uniform spacing. */ const float imf_log10_mass_bin_size = @@ -212,24 +211,21 @@ inline static void init_imf(struct feedback_props *restrict feedback_props) { (float)(feedback_props->n_imf_mass_bins - 1); /* Allocate IMF array */ - if (swift_memalign( - "imf-tables", (void **)&feedback_props->imf, - SWIFT_STRUCT_ALIGNMENT, - feedback_props->n_imf_mass_bins * sizeof(float)) != 0) + if (swift_memalign("imf-tables", (void **)&feedback_props->imf, + SWIFT_STRUCT_ALIGNMENT, + feedback_props->n_imf_mass_bins * sizeof(float)) != 0) error("Failed to allocate IMF bins table"); /* Allocate array to store IMF mass bins */ - if (swift_memalign( - "imf-tables", (void **)&feedback_props->imf_mass_bin, - SWIFT_STRUCT_ALIGNMENT, - feedback_props->n_imf_mass_bins * sizeof(float)) != 0) + if (swift_memalign("imf-tables", (void **)&feedback_props->imf_mass_bin, + SWIFT_STRUCT_ALIGNMENT, + feedback_props->n_imf_mass_bins * sizeof(float)) != 0) error("Failed to allocate IMF bins table"); /* Allocate array to store IMF mass bins in log10 space */ - if (swift_memalign( - "imf-tables", (void **)&feedback_props->imf_mass_bin_log10, - SWIFT_STRUCT_ALIGNMENT, - feedback_props->n_imf_mass_bins * sizeof(float)) != 0) + if (swift_memalign("imf-tables", (void **)&feedback_props->imf_mass_bin_log10, + SWIFT_STRUCT_ALIGNMENT, + feedback_props->n_imf_mass_bins * sizeof(float)) != 0) error("Failed to allocate IMF bins table"); /* Set IMF from Chabrier 2003 */ @@ -237,8 +233,7 @@ inline static void init_imf(struct feedback_props *restrict feedback_props) { for (int i = 0; i < feedback_props->n_imf_mass_bins; i++) { const float log10_mass_msun = - feedback_props->log10_imf_min_mass_msun + - i * imf_log10_mass_bin_size; + feedback_props->log10_imf_min_mass_msun + i * imf_log10_mass_bin_size; const float mass_msun = exp10f(log10_mass_msun); @@ -260,11 +255,10 @@ inline static void init_imf(struct feedback_props *restrict feedback_props) { } /* Normalize the IMF */ - const float norm = - integrate_imf(feedback_props->log10_imf_min_mass_msun, - feedback_props->log10_imf_max_mass_msun, - eagle_imf_integration_mass_weight, - /*(stellar_yields=)*/ NULL, feedback_props); + const float norm = integrate_imf(feedback_props->log10_imf_min_mass_msun, + feedback_props->log10_imf_max_mass_msun, + eagle_imf_integration_mass_weight, + /*(stellar_yields=)*/ NULL, feedback_props); for (int i = 0; i < feedback_props->n_imf_mass_bins; i++) feedback_props->imf[i] /= norm; @@ -281,8 +275,9 @@ inline static void init_imf(struct feedback_props *restrict feedback_props) { * @param star_properties the #stars_props data structure. * @return Mass of stars died up to that age in solar masses. */ -inline static float dying_mass_msun(const float age_Gyr, const float Z, - const struct feedback_props *feedback_props) { +inline static float dying_mass_msun( + const float age_Gyr, const float Z, + const struct feedback_props *feedback_props) { /* Pull out some common terms */ const double *lifetime_Z = feedback_props->lifetimes.metallicity; @@ -290,7 +285,7 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, double **const dying_times = feedback_props->lifetimes.dyingtime; const int n_Z = feedback_props->lifetimes.n_z; const int n_m = feedback_props->lifetimes.n_mass; - + /* Early abort? */ if (age_Gyr <= 0.f) { return feedback_props->imf_max_mass_msun; @@ -306,23 +301,22 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, /* Before start of the table */ Z_index = 0; Z_offset = 0.f; - + } else if (Z >= lifetime_Z[n_Z - 1]) { /* After end of the table */ Z_index = n_Z - 2; Z_offset = 1.f; - + } else { /* Normal case: Somewhere inside the table */ for (Z_index = 0; Z_index < n_Z - 1; Z_index++) { - if (lifetime_Z[Z_index + 1] > Z) - break; + if (lifetime_Z[Z_index + 1] > Z) break; } Z_offset = (Z - lifetime_Z[Z_index]) / - (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]); + (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]); } /* Check whether we are not beyond the table */ @@ -333,7 +327,7 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, /* Before start of the table */ time_index_lowZ = 0; time_offset_lowZ = 0.f; - + } else if (log10_age_yr <= dying_times[Z_index][n_m - 1]) { /* After end of the table */ @@ -349,7 +343,7 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, /* Before start of the table */ time_index_highZ = 0; time_offset_highZ = 0.f; - + } else if (log10_age_yr <= dying_times[Z_index + 1][n_m - 1]) { /* After end of the table */ @@ -361,9 +355,9 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, a solution for the two bounds */ int i = n_m; while (i >= 0 && (time_index_lowZ == -1 || time_index_highZ == -1)) { - + i--; - + if (dying_times[Z_index][i] >= log10_age_yr && time_index_lowZ == -1) { /* record index */ @@ -375,7 +369,7 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, (dying_times[Z_index][time_index_lowZ + 1] - dying_times[Z_index][time_index_lowZ]); } - + if (dying_times[Z_index + 1][i] >= log10_age_yr && time_index_highZ == -1) { /* record index */ @@ -383,15 +377,17 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, /* record distance from table element */ time_offset_highZ = - (log10_age_yr - dying_times[Z_index + 1][time_index_highZ]) / + (log10_age_yr - dying_times[Z_index + 1][time_index_highZ]) / (dying_times[Z_index + 1][time_index_highZ + 1] - dying_times[Z_index + 1][time_index_highZ]); } } /* And now interpolate the solution */ - const float mass_low_Z = interpolate_1d(lifetime_m, time_index_lowZ, time_offset_lowZ); - const float mass_high_Z = interpolate_1d(lifetime_m, time_index_highZ, time_offset_highZ); + const float mass_low_Z = + interpolate_1d(lifetime_m, time_index_lowZ, time_offset_lowZ); + const float mass_high_Z = + interpolate_1d(lifetime_m, time_index_highZ, time_offset_highZ); float mass = (1.f - Z_offset) * mass_low_Z + Z_offset * mass_high_Z; @@ -411,8 +407,9 @@ inline static float dying_mass_msun(const float age_Gyr, const float Z, * @param feedback_props the #feedback_props data structure. * @return The life time in Giga-years. */ -inline static float lifetime_in_Gyr(const float mass, const float Z, - const struct feedback_props *feedback_props) { +inline static float lifetime_in_Gyr( + const float mass, const float Z, + const struct feedback_props *feedback_props) { /* Pull out some common terms */ const double *lifetime_Z = feedback_props->lifetimes.metallicity; @@ -420,31 +417,30 @@ inline static float lifetime_in_Gyr(const float mass, const float Z, double **const dying_times = feedback_props->lifetimes.dyingtime; const int n_Z = feedback_props->lifetimes.n_z; const int n_m = feedback_props->lifetimes.n_mass; - + /* Calculate index along the mass axis */ int m_index; float m_offset; if (mass <= lifetime_m[0]) { - + /* Before start of the table */ m_index = 0; m_offset = 0.f; - + } else if (mass >= lifetime_m[n_m - 1]) { - + /* After end of the table */ m_index = n_m - 2; m_offset = 1.f; - + } else { /* Normal case: Somewhere inside the table */ for (m_index = 0; m_index < n_m - 1; m_index++) - if (lifetime_m[m_index + 1] > mass) - break; - + if (lifetime_m[m_index + 1] > mass) break; + m_offset = (mass - lifetime_m[m_index]) / - (lifetime_m[m_index + 1] - lifetime_m[m_index]); + (lifetime_m[m_index + 1] - lifetime_m[m_index]); } /* Calculate index along the metallicity axis */ @@ -455,30 +451,30 @@ inline static float lifetime_in_Gyr(const float mass, const float Z, /* Before start of the table */ Z_index = 0; Z_offset = 0.f; - + } else if (Z >= lifetime_Z[n_Z - 1]) { /* After end of the table */ Z_index = n_Z - 2; Z_offset = 1.f; - + } else { - + for (Z_index = 0; Z_index < n_Z - 1; Z_index++) - if (lifetime_Z[Z_index + 1] > Z) - break; + if (lifetime_Z[Z_index + 1] > Z) break; /* Normal case: Somewhere inside the table */ Z_offset = (Z - lifetime_Z[Z_index]) / - (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]); + (lifetime_Z[Z_index + 1] - lifetime_Z[Z_index]); } /* Interpolation of the table to get the time */ - const float log_time_years = interpolate_2d(dying_times, Z_index, m_index, Z_offset, m_offset); + const float log_time_years = + interpolate_2d(dying_times, Z_index, m_index, Z_offset, m_offset); /* Convert to Giga-years */ const float time_Gyr = exp10f(log_time_years - 9.f); - + return time_Gyr; } diff --git a/src/feedback/EAGLE/interpolate.h b/src/feedback/EAGLE/interpolate.h index a4b8ff1615a6146ef7084ea0f8c324d595a8c520..f9ef3f8fa23b6a60bb96cb4ea9e8f3decf5c8076 100644 --- a/src/feedback/EAGLE/interpolate.h +++ b/src/feedback/EAGLE/interpolate.h @@ -28,10 +28,8 @@ * @param i, j Indices of element of interest * @param Nx, Ny Sizes of array dimensions */ -__attribute__((always_inline)) static INLINE int row_major_index_2d(const int i, - const int j, - const int Nx, - const int Ny) { +__attribute__((always_inline)) static INLINE int row_major_index_2d( + const int i, const int j, const int Nx, const int Ny) { #ifdef SWIFT_DEBUG_CHECKS assert(i < Nx); assert(j < Ny); @@ -48,7 +46,8 @@ __attribute__((always_inline)) static INLINE int row_major_index_2d(const int i, * @param Nx, Ny, Nz Sizes of array dimensions */ __attribute__((always_inline)) static INLINE int row_major_index_3d( - const int i, const int j, const int k, const int Nx, const int Ny, const int Nz) { + const int i, const int j, const int k, const int Nx, const int Ny, + const int Nz) { #ifdef SWIFT_DEBUG_CHECKS assert(i < Nx); @@ -66,7 +65,8 @@ __attribute__((always_inline)) static INLINE int row_major_index_3d( * @param i index of cell to interpolate * @param dx offset within cell to interpolate */ -__attribute__((always_inline)) static INLINE double interpolate_1d(const double* table, const int i, const float dx) { +__attribute__((always_inline)) static INLINE double interpolate_1d( + const double* table, const int i, const float dx) { const float tx = 1.f - dx; @@ -82,8 +82,8 @@ __attribute__((always_inline)) static INLINE double interpolate_1d(const double* * @param dx row offset within cell to interpolate * @param dy column offset within cell to interpolate */ -__attribute__((always_inline)) static INLINE double interpolate_2d(double** table, const int i, const int j, const float dx, - const float dy) { +__attribute__((always_inline)) static INLINE double interpolate_2d( + double** table, const int i, const int j, const float dx, const float dy) { const float tx = 1.f - dx; const float ty = 1.f - dy; @@ -109,10 +109,11 @@ __attribute__((always_inline)) static INLINE double interpolate_2d(double** tabl * array_y to interpolate */ static INLINE double interpolate_1D_non_uniform(const double* array_x, - const double* array_y, const int size, + const double* array_y, + const int size, const double x) { #ifdef SWIFT_DEBUG_CHECKS - + /* Check that x within range of array_x */ if (x < array_x[0]) error("interpolating value less than array min. value %.5e array min %.5e", @@ -122,14 +123,14 @@ static INLINE double interpolate_1D_non_uniform(const double* array_x, "interpolating value greater than array max. value %.5e array max %.5e", x, array_x[size - 1]); #endif - + /* Find bin index and offset of x within array_x */ int index = 0; while (array_x[index] <= x) index++; const double offset = - (array_x[index] - x) / (array_x[index] - array_x[index - 1]); - + (array_x[index] - x) / (array_x[index] - array_x[index - 1]); + /* Interpolate array_y */ return offset * array_y[index - 1] + (1. - offset) * array_y[index]; } diff --git a/src/feedback/EAGLE/yield_tables.h b/src/feedback/EAGLE/yield_tables.h index cb84353fe3570ecf10bd64b52b229ce9abf49ea8..c26d0d38b48fcbdb7b79ac2610a3086cceb43ab8 100644 --- a/src/feedback/EAGLE/yield_tables.h +++ b/src/feedback/EAGLE/yield_tables.h @@ -183,14 +183,14 @@ inline static void read_yield_tables(struct feedback_props *feedback_props) { /* Flatten the temporary tables that were read, store in stars_props */ for (k = 0; k < feedback_props->SNII_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->SNII_n_z, - feedback_props->SNII_n_mass); + feedback_props->SNII_n_mass); feedback_props->yield_SNII.ejecta[flat_index] = temp_ejecta_SNII[k]; feedback_props->yield_SNII.total_metals[flat_index] = tempmet1[k]; for (j = 0; j < feedback_props->SNII_n_elements; j++) { - flat_index = row_major_index_3d( - i, j, k, feedback_props->SNII_n_z, feedback_props->SNII_n_elements, - feedback_props->SNII_n_mass); + flat_index = row_major_index_3d(i, j, k, feedback_props->SNII_n_z, + feedback_props->SNII_n_elements, + feedback_props->SNII_n_mass); feedback_props->yield_SNII.yield[flat_index] = temp_yield_SNII[j][k]; } } @@ -285,14 +285,14 @@ inline static void read_yield_tables(struct feedback_props *feedback_props) { /* Flatten the temporary tables that were read, store in stars_props */ for (k = 0; k < feedback_props->AGB_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->AGB_n_z, - feedback_props->AGB_n_mass); + feedback_props->AGB_n_mass); feedback_props->yield_AGB.ejecta[flat_index] = temp_ejecta_AGB[k]; feedback_props->yield_AGB.total_metals[flat_index] = tempmet2[k]; for (j = 0; j < feedback_props->AGB_n_elements; j++) { - flat_index = row_major_index_3d( - i, j, k, feedback_props->AGB_n_z, feedback_props->AGB_n_elements, - feedback_props->AGB_n_mass); + flat_index = row_major_index_3d(i, j, k, feedback_props->AGB_n_z, + feedback_props->AGB_n_elements, + feedback_props->AGB_n_mass); feedback_props->yield_AGB.yield[flat_index] = temp_yield_AGB[j][k]; } } @@ -347,7 +347,8 @@ inline static void read_yield_tables(struct feedback_props *feedback_props) { * * @param stars the #stars_props data struct to store the tables in */ -inline static void allocate_yield_tables(struct feedback_props *feedback_props) { +inline static void allocate_yield_tables( + struct feedback_props *feedback_props) { /* Allocate array to store SNIa yield tables */ if (swift_memalign("feedback-tables", (void **)&feedback_props->yields_SNIa, @@ -650,9 +651,9 @@ inline static void compute_yields(struct feedback_props *feedback_props) { feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]); } - flat_index_3d = row_major_index_3d( - i, elem, k, feedback_props->SNII_n_z, chemistry_element_count, - feedback_props->n_imf_mass_bins); + flat_index_3d = row_major_index_3d(i, elem, k, feedback_props->SNII_n_z, + chemistry_element_count, + feedback_props->n_imf_mass_bins); feedback_props->yield_SNII.yield_IMF_resampled[flat_index_3d] = exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result; } @@ -660,11 +661,11 @@ inline static void compute_yields(struct feedback_props *feedback_props) { for (int i = 0; i < feedback_props->SNII_n_z; i++) { for (int k = 0; k < feedback_props->n_imf_mass_bins; k++) { - flat_index_2d = row_major_index_2d( - i, k, feedback_props->SNII_n_z, feedback_props->n_imf_mass_bins); - flat_index_3d = row_major_index_3d( - i, elem, k, feedback_props->SNII_n_z, chemistry_element_count, - feedback_props->n_imf_mass_bins); + flat_index_2d = row_major_index_2d(i, k, feedback_props->SNII_n_z, + feedback_props->n_imf_mass_bins); + flat_index_3d = row_major_index_3d(i, elem, k, feedback_props->SNII_n_z, + chemistry_element_count, + feedback_props->n_imf_mass_bins); if (strcmp(chemistry_get_element_name(elem), "Hydrogen") != 0 || strcmp(chemistry_get_element_name(elem), "Helium") != 0) { feedback_props->yield_SNII @@ -737,7 +738,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { for (int i = 0; i < feedback_props->SNII_n_z; i++) { for (int k = 0; k < feedback_props->SNII_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->SNII_n_z, - feedback_props->SNII_n_mass); + feedback_props->SNII_n_mass); SNII_ejecta[k] = feedback_props->yield_SNII.ejecta[flat_index] * exp(M_LN10 * (-feedback_props->yield_SNII.mass[k])); } @@ -755,7 +756,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]); flat_index = row_major_index_2d(i, k, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); + feedback_props->n_imf_mass_bins); feedback_props->yield_SNII.ejecta_IMF_resampled[flat_index] = exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result; } @@ -765,7 +766,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { for (int i = 0; i < feedback_props->SNII_n_z; i++) { for (int k = 0; k < feedback_props->SNII_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->SNII_n_z, - feedback_props->SNII_n_mass); + feedback_props->SNII_n_mass); SNII_ejecta[k] = feedback_props->yield_SNII.total_metals[flat_index] * exp(M_LN10 * (-feedback_props->yield_SNII.mass[k])); } @@ -783,7 +784,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { feedback_props->SNII_n_mass, feedback_props->yield_mass_bins[k]); flat_index = row_major_index_2d(i, k, feedback_props->SNII_n_z, - feedback_props->n_imf_mass_bins); + feedback_props->n_imf_mass_bins); feedback_props->yield_SNII.total_metals_IMF_resampled[flat_index] = exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result; } @@ -793,7 +794,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { for (int i = 0; i < feedback_props->AGB_n_z; i++) { for (int k = 0; k < feedback_props->AGB_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->AGB_n_z, - feedback_props->AGB_n_mass); + feedback_props->AGB_n_mass); AGB_ejecta[k] = feedback_props->yield_AGB.ejecta[flat_index] / exp(M_LN10 * feedback_props->yield_AGB.mass[k]); } @@ -811,7 +812,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { feedback_props->AGB_n_mass, feedback_props->yield_mass_bins[k]); flat_index = row_major_index_2d(i, k, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); + feedback_props->n_imf_mass_bins); feedback_props->yield_AGB.ejecta_IMF_resampled[flat_index] = exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result; } @@ -820,7 +821,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { for (int i = 0; i < feedback_props->AGB_n_z; i++) { for (int k = 0; k < feedback_props->AGB_n_mass; k++) { flat_index = row_major_index_2d(i, k, feedback_props->AGB_n_z, - feedback_props->AGB_n_mass); + feedback_props->AGB_n_mass); AGB_ejecta[k] = feedback_props->yield_AGB.total_metals[flat_index] * exp(M_LN10 * (-feedback_props->yield_AGB.mass[k])); } @@ -838,7 +839,7 @@ inline static void compute_ejecta(struct feedback_props *feedback_props) { feedback_props->AGB_n_mass, feedback_props->yield_mass_bins[k]); flat_index = row_major_index_2d(i, k, feedback_props->AGB_n_z, - feedback_props->n_imf_mass_bins); + feedback_props->n_imf_mass_bins); feedback_props->yield_AGB.total_metals_IMF_resampled[flat_index] = exp(M_LN10 * feedback_props->yield_mass_bins[k]) * result; } diff --git a/src/runner_doiact_stars.h b/src/runner_doiact_stars.h index 3eb4a5526274824576c0983adb02afbe09018a9c..3dc1f1e5d07739ac1d8db6bc5bbbb3cb758b9835 100644 --- a/src/runner_doiact_stars.h +++ b/src/runner_doiact_stars.h @@ -99,7 +99,7 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) { /* Cosmological terms */ const float a = cosmo->a; const float H = cosmo->H; - + const int scount = c->stars.count; const int count = c->hydro.count; struct spart *restrict sparts = c->stars.parts; @@ -146,11 +146,11 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) { if (r2 < hig2) { IACT_STARS(r2, dx, hi, hj, si, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, cosmo, + feedback_props, xpj, ti_current); #endif } } /* loop over the parts in ci. */ @@ -187,7 +187,7 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci, /* Cosmological terms */ const float a = cosmo->a; const float H = cosmo->H; - + const int scount_i = ci->stars.count; const int count_j = cj->hydro.count; struct spart *restrict sparts_i = ci->stars.parts; @@ -242,15 +242,14 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci, if (r2 < hig2) { IACT_STARS(r2, dx, hi, hj, si, pj, a, H); - + #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, cosmo, + feedback_props, xpj, ti_current); #endif - } } /* loop over the parts in cj. */ } /* loop over the parts in ci. */ @@ -276,7 +275,7 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj, /* Cosmological terms */ const float a = cosmo->a; const float H = cosmo->H; - + /* Get the cutoff shift. */ double rshift = 0.0; for (int k = 0; k < 3; k++) rshift += shift[k] * runner_shift[sid][k]; @@ -409,13 +408,12 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj, IACT_STARS(r2, dx, hi, hj, spi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); #endif - } } /* loop over the parts in cj. */ } /* loop over the parts in ci. */ @@ -536,12 +534,12 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj, IACT_STARS(r2, dx, hj, hi, spj, pi, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hj, hi, spj, pi, cosmo, - feedback_props, xpi, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hj, hi, spj, pi, cosmo, + feedback_props, xpi, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hj, hi, spj, pi, cosmo, - feedback_props, xpi, ti_current); -#endif + runner_iact_nonsym_feedback_apply(r2, dx, hj, hi, spj, pi, cosmo, + feedback_props, xpi, ti_current); +#endif } } /* loop over the parts in ci. */ } /* loop over the parts in cj. */ @@ -635,7 +633,7 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, const double pjx = pj->x[0]; const double pjy = pj->x[1]; const double pjz = pj->x[2]; - const float hj = pj->h; + const float hj = pj->h; /* Compute the pairwise distance. */ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy), @@ -653,14 +651,14 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, /* Hit or miss? */ if (r2 < hig2) { IACT_STARS(r2, dx, hi, hj, spi, pj, a, H); - + #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); -#endif + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); +#endif } } /* loop over the parts in cj. */ } /* loop over the sparts in ci. */ @@ -695,8 +693,8 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, const double pjx = pj->x[0]; const double pjy = pj->x[1]; const double pjz = pj->x[2]; - const float hj = pj->h; - + const float hj = pj->h; + /* Compute the pairwise distance. */ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy), (float)(piz - pjz)}; @@ -715,12 +713,12 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, IACT_STARS(r2, dx, hi, hj, spi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); -#endif + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); +#endif } } /* loop over the parts in cj. */ } /* loop over the sparts in ci. */ @@ -762,7 +760,7 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci, const int count_j = cj->hydro.count; struct part *restrict parts_j = cj->hydro.parts; struct xpart *restrict xparts_j = cj->hydro.xparts; - + /* Early abort? */ if (count_j == 0) return; @@ -800,9 +798,9 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci, /* Compute the pairwise distance. */ float dx[3] = {(float)(pix - pjx), (float)(piy - pjy), - (float)(piz - pjz)}; + (float)(piz - pjz)}; const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; - + #ifdef SWIFT_DEBUG_CHECKS /* Check that particles have been drifted to the current time */ if (pj->ti_drift != e->ti_current) @@ -813,12 +811,12 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci, IACT_STARS(r2, dx, hi, hj, spi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, - feedback_props, xpj, ti_current); -#endif + runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, cosmo, + feedback_props, xpj, ti_current); +#endif } } /* loop over the parts in cj. */ } /* loop over the parts in ci. */ @@ -850,7 +848,7 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, /* Cosmological terms */ const float a = cosmo->a; const float H = cosmo->H; - + const int count_i = ci->hydro.count; struct part *restrict parts_j = ci->hydro.parts; struct xpart *restrict xparts_j = ci->hydro.xparts; @@ -901,13 +899,12 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci, if (r2 < hig2) { IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H); #if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY) - runner_iact_nonsym_feedback_density(r2, dx, hi, pj->h, spi, pj, cosmo, - feedback_props, xpj, ti_current); + runner_iact_nonsym_feedback_density(r2, dx, hi, pj->h, spi, pj, cosmo, + feedback_props, xpj, ti_current); #elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK) - runner_iact_nonsym_feedback_apply(r2, dx, hi, pj->h, spi, pj, cosmo, - feedback_props, xpj, ti_current); -#endif - + runner_iact_nonsym_feedback_apply(r2, dx, hi, pj->h, spi, pj, cosmo, + feedback_props, xpj, ti_current); +#endif } } /* loop over the parts in cj. */ } /* loop over the parts in ci. */ diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h index 38a64d22553caa370dd96abefea0ba772d423330..9d20fd5bb535acb25d1fb24bbef61d88cc5a0b93 100644 --- a/src/stars/Default/stars_iact.h +++ b/src/stars/Default/stars_iact.h @@ -33,9 +33,10 @@ * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_stars_density( - float r2, const float *dx, float hi, float hj, struct spart *restrict si, - const struct part *restrict pj, const float a, const float H) { +runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj, + struct spart *restrict si, + const struct part *restrict pj, const float a, + const float H) { float wi, wi_dx; @@ -75,10 +76,10 @@ runner_iact_nonsym_stars_density( * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_stars_feedback( - float r2, const float *dx, float hi, float hj, - const struct spart *restrict si, struct part *restrict pj, - const float a, const float H) { +runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj, + const struct spart *restrict si, + struct part *restrict pj, const float a, + const float H) { const float mj = hydro_get_mass(pj); const float rhoj = hydro_get_comoving_density(pj); diff --git a/src/stars/EAGLE/stars_iact.h b/src/stars/EAGLE/stars_iact.h index 5c82e7777a9103d7c042405476d83c1f84523f62..7958eaec563778221dc8f7a69252890f42dd9a6d 100644 --- a/src/stars/EAGLE/stars_iact.h +++ b/src/stars/EAGLE/stars_iact.h @@ -34,9 +34,10 @@ * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_stars_density( - float r2, const float *dx, float hi, float hj, struct spart *restrict si, - const struct part *restrict pj, const float a, const float H) { +runner_iact_nonsym_stars_density(float r2, const float *dx, float hi, float hj, + struct spart *restrict si, + const struct part *restrict pj, const float a, + const float H) { float wi, wi_dx; @@ -81,10 +82,10 @@ runner_iact_nonsym_stars_density( * @param H Current Hubble parameter. */ __attribute__((always_inline)) INLINE static void -runner_iact_nonsym_stars_feedback( - float r2, const float *dx, float hi, float hj, - const struct spart *restrict si, struct part *restrict pj, - const float a, const float H) { +runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj, + const struct spart *restrict si, + struct part *restrict pj, const float a, + const float H) { #ifdef DEBUG_INTERACTIONS_STARS /* Update ngb counters */ @@ -92,7 +93,6 @@ runner_iact_nonsym_stars_feedback( si->ids_ngbs_feedback[si->num_ngb_feedback] = pj->id; ++si->num_ngb_feedback; #endif - } #endif /* SWIFT_EAGLE_STARS_IACT_H */ diff --git a/src/stars/EAGLE/stars_part.h b/src/stars/EAGLE/stars_part.h index 7bf82d03ecd056e08a339fcca59f891ca74b7f80..047f1d860476ed0c55dbc2d61ffc177e098bb0c3 100644 --- a/src/stars/EAGLE/stars_part.h +++ b/src/stars/EAGLE/stars_part.h @@ -100,7 +100,7 @@ struct spart { /*! Feedback structure */ struct feedback_part_data feedback_data; - + /*! Tracer structure */ struct tracers_xpart_data tracers_data; diff --git a/src/tools.c b/src/tools.c index 591bcdccc22ba001760ca4d5f078aea9e5fe6246..b4dc3a1ca2d27c6eacbd8cb55c95614468257e19 100644 --- a/src/tools.c +++ b/src/tools.c @@ -644,7 +644,7 @@ void self_all_force(struct runner *r, struct cell *ci) { } void self_all_stars_density(struct runner *r, struct cell *ci) { - + float r2, hi, hj, hig2, dxi[3]; struct spart *spi; struct part *pj;