Skip to content
Snippets Groups Projects
Commit 07fe0a7c authored by Zachary Rey's avatar Zachary Rey
Browse files

uncomment

parent cfda4b97
No related branches found
No related tags found
No related merge requests found
...@@ -102,7 +102,7 @@ void feedback_update_part(struct part* p, struct xpart* xp, ...@@ -102,7 +102,7 @@ void feedback_update_part(struct part* p, struct xpart* xp,
hydro_get_physical_internal_energy(p, xp, cosmo) * old_mass / new_mass; hydro_get_physical_internal_energy(p, xp, cosmo) * old_mass / new_mass;
const float u_new = u + xp->feedback_data.delta_u; const float u_new = u + xp->feedback_data.delta_u;
message("The particle %lld has a new energy : %lf compared to the old energy : %lf", p->id, u_new, u); //message("The particle %lld has a new energy : %lf compared to the old energy : %lf", p->id, u_new, u);
hydro_set_physical_internal_energy(p, xp, cosmo, u_new); hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, u_new); hydro_set_drifted_physical_internal_energy(p, cosmo, pressure_floor, u_new);
...@@ -272,7 +272,7 @@ void feedback_will_do_feedback( ...@@ -272,7 +272,7 @@ void feedback_will_do_feedback(
/* TODO: Do we want to multiply pre-SN energy by the efficiency? */ /* TODO: Do we want to multiply pre-SN energy by the efficiency? */
sp->feedback_data.preSN.energy_ejected *= feedback_props->preSN_efficiency; sp->feedback_data.preSN.energy_ejected *= feedback_props->preSN_efficiency;
message("Energy of preSN ejected with efficiency coefficient = %e (in internal units)", sp->feedback_data.preSN.energy_ejected); // message("Energy of preSN ejected with efficiency coefficient = %e (in internal units)", sp->feedback_data.preSN.energy_ejected);
/* TODO: See if we need to add something about pre-SN */ /* TODO: See if we need to add something about pre-SN */
/* Set the particle as doing some feedback */ /* Set the particle as doing some feedback */
......
...@@ -143,7 +143,7 @@ runner_iact_nonsym_feedback_apply( ...@@ -143,7 +143,7 @@ runner_iact_nonsym_feedback_apply(
if (e_preSN != 0.0) { if (e_preSN != 0.0) {
/* Energy received */ /* Energy received */
const double du = (e_preSN) * weight / new_mass; const double du = (e_preSN) * weight / new_mass;
message("the received energy of particle %lld is : %e",pj->id,du); //message("the received energy of particle %lld is : %e",pj->id,du);
xpj->feedback_data.delta_u += du; xpj->feedback_data.delta_u += du;
} }
......
...@@ -363,13 +363,13 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -363,13 +363,13 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
log_metallicity = log10(metallicity / 0.02); log_metallicity = log10(metallicity / 0.02);
} }
message("The log metallicity in solar metallicity is : %f", log_metallicity); //message("The log metallicity in solar metallicity is : %f", log_metallicity);
/* If the star particle the calculation is straight forward */ /* If the star particle the calculation is straight forward */
if (sp->star_type == single_star) { if (sp->star_type == single_star) {
/* If single star, only the mass of the star to consider */ /* If single star, only the mass of the star to consider */
const float log_m = log10(m_beg_lim); const float log_m = log10(m_beg_lim);
message("The solar mass considered is %g", m_beg_lim); //message("The solar mass considered is %g", m_beg_lim);
/* Stellar winds */ /* Stellar winds */
/* Compute the mass-loss */ /* Compute the mass-loss */
sp->feedback_data.preSN.mass_loss = stellar_wind_get_ejected_mass(log_metallicity, log_m); sp->feedback_data.preSN.mass_loss = stellar_wind_get_ejected_mass(log_metallicity, log_m);
...@@ -380,10 +380,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -380,10 +380,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
/* Stellar winds contribution */ /* Stellar winds contribution */
sp->feedback_data.preSN.energy_ejected = stellar_wind_get_energy_dot(sp->feedback_data.preSN.mass_loss,v_infinity); sp->feedback_data.preSN.energy_ejected = stellar_wind_get_energy_dot(sp->feedback_data.preSN.mass_loss,v_infinity);
message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e", // message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e",
sp->feedback_data.preSN.mass_loss, // sp->feedback_data.preSN.mass_loss,
v_infinity, // v_infinity,
sp->feedback_data.preSN.energy_ejected); // sp->feedback_data.preSN.energy_ejected);
/* Radiation ? */ /* Radiation ? */
...@@ -410,7 +410,7 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -410,7 +410,7 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
float mass_loss =0; float mass_loss =0;
float v_infty = 0; float v_infty = 0;
message("The imf star has an imf mass of %e; The m_end_lim is %e;", imf_m, m_end_lim); //message("The imf star has an imf mass of %e; The m_end_lim is %e;", imf_m, m_end_lim);
/* Calculate the part of the imf which is not considered yet for supernovae*/ /* Calculate the part of the imf which is not considered yet for supernovae*/
while (imf_m + dM < m_end_lim){ while (imf_m + dM < m_end_lim){
...@@ -423,10 +423,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -423,10 +423,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m; energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m;
imf_m = imf_m + dM; imf_m = imf_m + dM;
dM = imf_m / 10; dM = imf_m / 10;
message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e", // message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e",
mass_loss, // mass_loss,
v_infty, // v_infty,
energy_dot); // energy_dot);
} }
if (imf_m < m_end_lim){ if (imf_m < m_end_lim){
log_m = log10((imf_m + m_end_lim) / 2); log_m = log10((imf_m + m_end_lim) / 2);
...@@ -436,10 +436,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -436,10 +436,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
v_infty = stellar_wind_get_wind_velocity(log_metallicity, log_m); v_infty = stellar_wind_get_wind_velocity(log_metallicity, log_m);
sp->feedback_data.preSN.mass_loss += mass_loss * N_star_m; sp->feedback_data.preSN.mass_loss += mass_loss * N_star_m;
energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m; energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m;
message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e", // message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e",
mass_loss, // mass_loss,
v_infty, // v_infty,
energy_dot); // energy_dot);
} }
/* Then consider the imf part which could provoke supernovae*/ /* Then consider the imf part which could provoke supernovae*/
...@@ -455,10 +455,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -455,10 +455,10 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m; energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m;
imf_m = imf_m + dM; imf_m = imf_m + dM;
dM = imf_m / 10; dM = imf_m / 10;
message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e", // message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e",
mass_loss, // mass_loss,
v_infty, // v_infty,
energy_dot); // energy_dot);
} }
if (imf_m < m_beg_lim){ if (imf_m < m_beg_lim){
log_m = log10((imf_m + m_beg_lim) / 2); log_m = log10((imf_m + m_beg_lim) / 2);
...@@ -468,14 +468,14 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const ...@@ -468,14 +468,14 @@ void stellar_evolution_compute_preSN_properties(struct spart* restrict sp, const
v_infty = stellar_wind_get_wind_velocity(log_metallicity, log_m); v_infty = stellar_wind_get_wind_velocity(log_metallicity, log_m);
sp->feedback_data.preSN.mass_loss += mass_loss * N_star_m; sp->feedback_data.preSN.mass_loss += mass_loss * N_star_m;
energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m; energy_dot += stellar_wind_get_energy_dot(mass_loss, v_infty) * N_star_m;
message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e", // message("The Mass-Loss = %e; The wind velocity = %e; The energy_dot = %e",
mass_loss, // mass_loss,
v_infty, // v_infty,
energy_dot); // energy_dot);
} }
sp->feedback_data.preSN.energy_ejected = energy_dot; sp->feedback_data.preSN.energy_ejected = energy_dot;
message("The energy ejected of imf star = %e (in erg)", energy_dot); //message("The energy ejected of imf star = %e (in erg)", energy_dot);
} }
} }
...@@ -587,8 +587,6 @@ void stellar_evolution_evolve_spart( ...@@ -587,8 +587,6 @@ void stellar_evolution_evolve_spart(
return; return;
} }
message("Is this even called ???");
/* TODO: Pre-SN feedback */ /* TODO: Pre-SN feedback */
/* Note: You can update the function parameters as needed. */ /* Note: You can update the function parameters as needed. */
stellar_evolution_compute_preSN_feedback_spart(sp, sm, cosmo, us, phys_const, ti_begin, stellar_evolution_compute_preSN_feedback_spart(sp, sm, cosmo, us, phys_const, ti_begin,
...@@ -1196,7 +1194,7 @@ void stellar_evolution_compute_preSN_feedback_individual_star(struct spart* rest ...@@ -1196,7 +1194,7 @@ void stellar_evolution_compute_preSN_feedback_individual_star(struct spart* rest
const double star_age_beg_step, const double dt) { const double star_age_beg_step, const double dt) {
/* TODO */ /* TODO */
/* TODO erase this debug line*/ /* TODO erase this debug line*/
message("Computing individual preSN feedback for stellar particle : %lld", sp->id); //message("Computing individual preSN feedback for stellar particle : %lld", sp->id);
/* Check that this function is called for individual stars (REDUNDANT) */ /* Check that this function is called for individual stars (REDUNDANT) */
if (sp->star_type != single_star) { if (sp->star_type != single_star) {
...@@ -1241,13 +1239,13 @@ void stellar_evolution_compute_preSN_feedback_individual_star(struct spart* rest ...@@ -1241,13 +1239,13 @@ void stellar_evolution_compute_preSN_feedback_individual_star(struct spart* rest
sp->feedback_data.preSN.energy_ejected *= feedback_duration_yr; sp->feedback_data.preSN.energy_ejected *= feedback_duration_yr;
message("The energy amount to eject is : %lf (in ergs)", sp->feedback_data.preSN.energy_ejected); //message("The energy amount to eject is : %lf (in ergs)", sp->feedback_data.preSN.energy_ejected);
/* convert to internal units */ /* convert to internal units */
sp->feedback_data.preSN.mass_loss *= phys_const->const_solar_mass; sp->feedback_data.preSN.mass_loss *= phys_const->const_solar_mass;
sp->feedback_data.preSN.energy_ejected /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY); sp->feedback_data.preSN.energy_ejected /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
message("The conversion factor for energy : %e (in internal units)", units_cgs_conversion_factor(us, UNIT_CONV_ENERGY)); //message("The conversion factor for energy : %e (in internal units)", units_cgs_conversion_factor(us, UNIT_CONV_ENERGY));
message("The energy amount to eject is : %e (in internal units)", sp->feedback_data.preSN.energy_ejected); //message("The energy amount to eject is : %e (in internal units)", sp->feedback_data.preSN.energy_ejected);
/* maybe we want to consider also the radiation contribution */ /* maybe we want to consider also the radiation contribution */
...@@ -1311,7 +1309,7 @@ void stellar_evolution_compute_preSN_feedback_spart( ...@@ -1311,7 +1309,7 @@ void stellar_evolution_compute_preSN_feedback_spart(
m_end_step = max(m_end_step, sm->imf.mass_min); m_end_step = max(m_end_step, sm->imf.mass_min);
m_beg_step = min(m_beg_step, sm->imf.mass_max); m_beg_step = min(m_beg_step, sm->imf.mass_max);
message("The m_end_step is : %e and the m_beg_step is : %e (in ? units)",m_end_step,m_beg_step); // message("The m_end_step is : %e and the m_beg_step is : %e (in ? units)",m_end_step,m_beg_step);
/* considering only the "alive" part of the IMF, i.e., we stop only if we are currently below the IMF */ /* considering only the "alive" part of the IMF, i.e., we stop only if we are currently below the IMF */
if (m_beg_step < sm->imf.mass_min) return; if (m_beg_step < sm->imf.mass_min) return;
...@@ -1326,7 +1324,7 @@ void stellar_evolution_compute_preSN_feedback_spart( ...@@ -1326,7 +1324,7 @@ void stellar_evolution_compute_preSN_feedback_spart(
if (sp->star_type == star_population_continuous_IMF) { if (sp->star_type == star_population_continuous_IMF) {
/* If it's not time yet for feedback, exit. Notice that both masses are in /* If it's not time yet for feedback, exit. Notice that both masses are in
solar mass. */ solar mass. */
message("The minimal imf discret mass is %e (in Msun)",sm->imf.minimal_discrete_mass_Msun); //message("The minimal imf discret mass is %e (in Msun)",sm->imf.minimal_discrete_mass_Msun);
/* If we are in a case where /* If we are in a case where
m_beg_step and/or m_end_step > minimal_discrete_mass_Msun, m_beg_step and/or m_end_step > minimal_discrete_mass_Msun,
...@@ -1349,13 +1347,13 @@ void stellar_evolution_compute_preSN_feedback_spart( ...@@ -1349,13 +1347,13 @@ void stellar_evolution_compute_preSN_feedback_spart(
stellar_evolution_compute_preSN_properties(sp, sm, phys_const,m_beg_step, m_end_step, m_init); stellar_evolution_compute_preSN_properties(sp, sm, phys_const,m_beg_step, m_end_step, m_init);
sp->feedback_data.preSN.energy_ejected *= dt / phys_const->const_year; sp->feedback_data.preSN.energy_ejected *= dt / phys_const->const_year;
message("The energy amount to eject is : %e (in ergs)", sp->feedback_data.preSN.energy_ejected); //message("The energy amount to eject is : %e (in ergs)", sp->feedback_data.preSN.energy_ejected);
/* convert to internal units */ /* convert to internal units */
sp->feedback_data.preSN.mass_loss *= phys_const->const_solar_mass; sp->feedback_data.preSN.mass_loss *= phys_const->const_solar_mass;
sp->feedback_data.preSN.energy_ejected /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY); sp->feedback_data.preSN.energy_ejected /= units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
message("The conversion factor for energy : %e (in internal units)", units_cgs_conversion_factor(us, UNIT_CONV_ENERGY)); //message("The conversion factor for energy : %e (in internal units)", units_cgs_conversion_factor(us, UNIT_CONV_ENERGY));
message("The energy amount to eject is : %e (in internal units)", sp->feedback_data.preSN.energy_ejected); //message("The energy amount to eject is : %e (in internal units)", sp->feedback_data.preSN.energy_ejected);
//TODO the same thought than for individual stars //TODO the same thought than for individual stars
} }
\ No newline at end of file
...@@ -97,7 +97,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){ ...@@ -97,7 +97,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){
/* If the star is lower than a limit mass, calculate the function*/ /* If the star is lower than a limit mass, calculate the function*/
if(log_m < stellar_wind_x0){ if(log_m < stellar_wind_x0){
for (int i=0; i < 4; i++){ for (int i=0; i < 4; i++){
message("Inside wind velocity, calculate parameter %i is = %f", i, calculate_b_parameter(log_Z,coeffs[i])); //message("Inside wind velocity, calculate parameter %i is = %f", i, calculate_b_parameter(log_Z,coeffs[i]));
wind_velocity += calculate_b_parameter(log_Z,coeffs[i]) * pow(log_m,i); wind_velocity += calculate_b_parameter(log_Z,coeffs[i]) * pow(log_m,i);
} }
/* Else, the function will be only a linear relation between the initial mass and the function calculated for the limit mass */ /* Else, the function will be only a linear relation between the initial mass and the function calculated for the limit mass */
...@@ -107,7 +107,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){ ...@@ -107,7 +107,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){
for (int i=0; i < 4; i++){ for (int i=0; i < 4; i++){
A0 += calculate_b_parameter(log_Z,coeffs[i]) * pow(stellar_wind_x0,i); A0 += calculate_b_parameter(log_Z,coeffs[i]) * pow(stellar_wind_x0,i);
message("Inside wind velocity, calculate parameter %i is = %f", i, calculate_b_parameter(log_Z,coeffs[i])); //message("Inside wind velocity, calculate parameter %i is = %f", i, calculate_b_parameter(log_Z,coeffs[i]));
// The derivative take one less step in the loop // The derivative take one less step in the loop
if (i == 0){ if (i == 0){
continue; continue;
...@@ -118,7 +118,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){ ...@@ -118,7 +118,7 @@ float stellar_wind_get_wind_velocity(const float log_Z,const float log_m){
wind_velocity += A0 + dLogA0 * (log_m - stellar_wind_x0); wind_velocity += A0 + dLogA0 * (log_m - stellar_wind_x0);
} }
message("Inside wind velocity computation : V_infty = %f", wind_velocity); //message("Inside wind velocity computation : V_infty = %f", wind_velocity);
return exp10(wind_velocity); return exp10(wind_velocity);
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment