Commit 1d676231 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

GEAR: cleanup

parent d5f2dbb1
...@@ -21,7 +21,6 @@ ...@@ -21,7 +21,6 @@
/** /**
* @file src/chemistry/GEAR/chemistry.h * @file src/chemistry/GEAR/chemistry.h
* @brief Empty infrastructure for the cases without chemistry function
*/ */
/* Some standard headers. */ /* Some standard headers. */
......
...@@ -655,8 +655,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const, ...@@ -655,8 +655,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
const struct part* restrict p, const struct part* restrict p,
struct xpart* restrict xp) { struct xpart* restrict xp) {
error("TODO: use energy after adiabatic cooling");
/* set current time */ /* set current time */
code_units units = cooling->units; code_units units = cooling->units;
...@@ -711,31 +709,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const, ...@@ -711,31 +709,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
return cooling_time; return cooling_time;
} }
/**
* @brief Apply the cooling to a particle.
*
* Depending on the task order, you may wish to either
* cool down the particle immediately or do it during the drift.
*
* @param p Pointer to the particle data.
* @param xp Pointer to the xparticle data.
* @param cosmo The current cosmological model.
* @param cooling_du_dt Time derivative of the cooling.
* @param u_new Internal energy after the cooling.
*/
void cooling_apply(struct part* restrict p, struct xpart* restrict xp,
const struct cosmology* restrict cosmo, float cooling_du_dt,
gr_float u_new) {
#ifdef TASK_ORDER_GEAR
/* Cannot use du / dt as it will be erased before being used */
hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
hydro_set_drifted_physical_internal_energy(p, cosmo, u_new);
#else
hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
#endif
}
/** /**
* @brief Apply the cooling function to a particle. * @brief Apply the cooling function to a particle.
* *
...@@ -797,13 +770,6 @@ void cooling_cool_part(const struct phys_const* restrict phys_const, ...@@ -797,13 +770,6 @@ void cooling_cool_part(const struct phys_const* restrict phys_const,
/* Calculate the cooling rate */ /* Calculate the cooling rate */
float cool_du_dt = (u_new - u_ad_before) / dt_therm; float cool_du_dt = (u_new - u_ad_before) / dt_therm;
#ifdef TASK_ORDER_GEAR
/* Set the energy */
hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
hydro_set_drifted_physical_internal_energy(p, cosmo, u_new);
#endif
/* Check that the energy stays above the limits if the time step increase by 2 */ /* Check that the energy stays above the limits if the time step increase by 2 */
/* Hydro */ /* Hydro */
double u_ad = u_new + hydro_du_dt * dt_therm; double u_ad = u_new + hydro_du_dt * dt_therm;
......
...@@ -105,9 +105,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const, ...@@ -105,9 +105,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
const struct cosmology* restrict cosmo, const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling, const struct cooling_function_data* restrict cooling,
const struct part* restrict p, struct xpart* restrict xp); const struct part* restrict p, struct xpart* restrict xp);
void cooling_apply(struct part* restrict p, struct xpart* restrict xp,
const struct cosmology* restrict cosmo, float cooling_du_dt,
gr_float u_new);
void cooling_cool_part(const struct phys_const* restrict phys_const, void cooling_cool_part(const struct phys_const* restrict phys_const,
const struct unit_system* restrict us, const struct unit_system* restrict us,
......
...@@ -115,13 +115,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx, ...@@ -115,13 +115,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
/* Mass received */ /* Mass received */
const double m_ej = si->feedback_data.mass_ejected; const double m_ej = si->feedback_data.mass_ejected;
// TODO compute inverse before feedback loop
const double weight = mj * wi * si_inv_weight; const double weight = mj * wi * si_inv_weight;
const double dm = m_ej * weight; const double dm = m_ej * weight;
const double new_mass = mj + dm; const double new_mass = mj + dm;
/* Energy received */ /* Energy received */
// TODO compute inverse before feedback loop
const double du = e_sn * weight / new_mass; const double du = e_sn * weight / new_mass;
xpj->feedback_data.delta_mass += dm; xpj->feedback_data.delta_mass += dm;
......
...@@ -63,8 +63,6 @@ __attribute__((always_inline)) INLINE static void feedback_props_print( ...@@ -63,8 +63,6 @@ __attribute__((always_inline)) INLINE static void feedback_props_print(
/** /**
* @brief Initialize the global properties of the feedback scheme. * @brief Initialize the global properties of the feedback scheme.
* *
* By default, takes the values provided by the hydro.
*
* @param fp The #feedback_props. * @param fp The #feedback_props.
* @param phys_const The physical constants in the internal unit system. * @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system. * @param us The internal unit system.
......
...@@ -19,6 +19,9 @@ ...@@ -19,6 +19,9 @@
#ifndef SWIFT_GEAR_INTERPOLATION_H #ifndef SWIFT_GEAR_INTERPOLATION_H
#define SWIFT_GEAR_INTERPOLATION_H #define SWIFT_GEAR_INTERPOLATION_H
/**
* @brief Type of boundary condition available.
*/
enum interpolate_boundary_condition { enum interpolate_boundary_condition {
/* No extrapolation => raise errors */ /* No extrapolation => raise errors */
boundary_condition_error, boundary_condition_error,
...@@ -33,6 +36,9 @@ enum interpolate_boundary_condition { ...@@ -33,6 +36,9 @@ enum interpolate_boundary_condition {
boundary_condition_const, boundary_condition_const,
}; };
/**
* @brief Structure for the interpolation.
*/
struct interpolation_1d { struct interpolation_1d {
/* Data to interpolate */ /* Data to interpolate */
float *data; float *data;
...@@ -53,8 +59,6 @@ struct interpolation_1d { ...@@ -53,8 +59,6 @@ struct interpolation_1d {
/** /**
* @brief Initialize the #interpolation_1d. * @brief Initialize the #interpolation_1d.
* *
* Assumes x are linear in log.
*
* @params interp The #interpolation_1d. * @params interp The #interpolation_1d.
* @params xmin Minimal value of x (in log). * @params xmin Minimal value of x (in log).
* @params xmax Maximal value of x (in log). * @params xmax Maximal value of x (in log).
......
...@@ -111,7 +111,7 @@ void stellar_evolution_compute_continuous_feedback_properties( ...@@ -111,7 +111,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
/* SNII */ /* SNII */
const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed_from_integral( const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
&sm->snii, log_m_end_step, log_m_beg_step); &sm->snii, log_m_end_step, log_m_beg_step);
/* Sum the contributions from SNIa and SNII */ /* Sum the contributions from SNIa and SNII */
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
...@@ -136,7 +136,7 @@ void stellar_evolution_compute_continuous_feedback_properties( ...@@ -136,7 +136,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
snii_yields); snii_yields);
/* Compute the mass fraction of non processed elements */ /* Compute the mass fraction of non processed elements */
const float non_processed = supernovae_ii_get_ejected_mass_fraction_from_integral( const float non_processed = supernovae_ii_get_ejected_mass_fraction_non_processed_from_integral(
&sm->snii, log_m_end_step, log_m_beg_step); &sm->snii, log_m_end_step, log_m_beg_step);
/* Set the yields */ /* Set the yields */
...@@ -222,7 +222,7 @@ void stellar_evolution_compute_discrete_feedback_properties( ...@@ -222,7 +222,7 @@ void stellar_evolution_compute_discrete_feedback_properties(
/* Compute the mass fraction of non processed elements */ /* Compute the mass fraction of non processed elements */
const float non_processed = const float non_processed =
supernovae_ii_get_ejected_mass_fraction_from_raw( supernovae_ii_get_ejected_mass_fraction_non_processed_from_raw(
&sm->snii, log_m_avg); &sm->snii, log_m_avg);
/* Set the yields */ /* Set the yields */
......
...@@ -117,12 +117,12 @@ struct supernovae_ii { ...@@ -117,12 +117,12 @@ struct supernovae_ii {
struct { struct {
/*! Mass fraction of metals ejected by a supernovae. */ /*! Mass fraction of metals ejected by a supernovae. */
struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT]; struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
/*! Total mass fraction ejected. */ /*! Total mass fraction ejected. */
struct interpolation_1d ejected_mass_processed; struct interpolation_1d ejected_mass_processed;
/*! Mass fraction ejected and not processed (=> with the star metallicity). */ /*! Mass fraction ejected and not processed (=> with the star metallicity). */
struct interpolation_1d ejected_mass; struct interpolation_1d ejected_mass_non_processed;
} raw; } raw;
/*! Yields integrated */ /*! Yields integrated */
...@@ -130,12 +130,12 @@ struct supernovae_ii { ...@@ -130,12 +130,12 @@ struct supernovae_ii {
/*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
*/ */
struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT]; struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
/*! Total mass fraction ejected (integrated over the IMF) */ /*! Total mass fraction ejected (integrated over the IMF) */
struct interpolation_1d ejected_mass_processed; struct interpolation_1d ejected_mass_processed;
/*! Mass fraction ejected and not processed (=> with the star metallicity) */ /*! Mass fraction ejected and not processed (=> with the star metallicity) */
struct interpolation_1d ejected_mass; struct interpolation_1d ejected_mass_non_processed;
} integrated; } integrated;
/*! Minimal mass for a SNII */ /*! Minimal mass for a SNII */
......
...@@ -123,36 +123,36 @@ void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float l ...@@ -123,36 +123,36 @@ void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float l
}; };
/** /**
* @brief Get the ejected mass per mass unit. * @brief Get the ejected mass (non processed) per mass unit.
* *
* @param snii The #supernovae_ii model. * @param snii The #supernovae_ii model.
* @param m1 The lower mass in log. * @param m1 The lower mass in log.
* @param m2 The upper mass in log. * @param m2 The upper mass in log.
* *
* @return mass_ejected_processed The mass of processsed elements. * @return mass_ejected_processed The mass of non processsed elements.
*/ */
float supernovae_ii_get_ejected_mass_fraction_from_integral(const struct supernovae_ii *snii, float supernovae_ii_get_ejected_mass_fraction_non_processed_from_integral(const struct supernovae_ii *snii,
float log_m1, float log_m2) { float log_m1, float log_m2) {
float mass_ejected_1 = interpolate_1d(&snii->integrated.ejected_mass, log_m1); float mass_ejected_1 = interpolate_1d(&snii->integrated.ejected_mass_non_processed, log_m1);
float mass_ejected_2 = interpolate_1d(&snii->integrated.ejected_mass, log_m2); float mass_ejected_2 = interpolate_1d(&snii->integrated.ejected_mass_non_processed, log_m2);
return mass_ejected_2 - mass_ejected_1; return mass_ejected_2 - mass_ejected_1;
}; };
/** /**
* @brief Get the ejected mass per mass unit. * @brief Get the ejected mass (non processed) per mass unit.
* *
* @param snii The #supernovae_ii model. * @param snii The #supernovae_ii model.
* @param log_m The mass in log. * @param log_m The mass in log.
* *
* @return The mass of processsed elements. * @return The mass of non processsed elements.
*/ */
float supernovae_ii_get_ejected_mass_fraction_from_raw(const struct supernovae_ii *snii, float supernovae_ii_get_ejected_mass_fraction_non_processed_from_raw(const struct supernovae_ii *snii,
float log_m) { float log_m) {
return interpolate_1d(&snii->raw.ejected_mass, log_m); return interpolate_1d(&snii->raw.ejected_mass_non_processed, log_m);
}; };
/** /**
...@@ -307,8 +307,8 @@ void supernovae_ii_read_yields(struct supernovae_ii *snii, ...@@ -307,8 +307,8 @@ void supernovae_ii_read_yields(struct supernovae_ii *snii,
"Ej", &previous_count, interpolation_size); "Ej", &previous_count, interpolation_size);
/* Read the mass ejected of non processed gas */ /* Read the mass ejected of non processed gas */
supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass, supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass_non_processed,
&snii->integrated.ejected_mass, &snii->integrated.ejected_mass_non_processed,
phys_const, sm, group_id, "Ejnp", phys_const, sm, group_id, "Ejnp",
&previous_count, interpolation_size); &previous_count, interpolation_size);
...@@ -449,16 +449,16 @@ void supernovae_ii_dump(const struct supernovae_ii *snii, FILE *stream, ...@@ -449,16 +449,16 @@ void supernovae_ii_dump(const struct supernovae_ii *snii, FILE *stream,
} }
/*! Dump the non processed mass (integrated). */ /*! Dump the non processed mass (integrated). */
if (snii->integrated.ejected_mass.data != NULL) { if (snii->integrated.ejected_mass_non_processed.data != NULL) {
restart_write_blocks((void *)snii->integrated.ejected_mass.data, restart_write_blocks((void *)snii->integrated.ejected_mass_non_processed.data,
sizeof(float), snii->integrated.ejected_mass.N, stream, sizeof(float), snii->integrated.ejected_mass_non_processed.N, stream,
"non_processed_mass_int", "non_processed_mass_int"); "non_processed_mass_int", "non_processed_mass_int");
} }
/*! Dump the non processed mass (raw). */ /*! Dump the non processed mass (raw). */
if (snii->raw.ejected_mass.data != NULL) { if (snii->raw.ejected_mass_non_processed.data != NULL) {
restart_write_blocks((void *)snii->raw.ejected_mass.data, restart_write_blocks((void *)snii->raw.ejected_mass_non_processed.data,
sizeof(float), snii->raw.ejected_mass.N, stream, sizeof(float), snii->raw.ejected_mass_non_processed.N, stream,
"non_processed_mass_raw", "non_processed_mass_raw"); "non_processed_mass_raw", "non_processed_mass_raw");
} }
...@@ -550,28 +550,28 @@ void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream, ...@@ -550,28 +550,28 @@ void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream,
} }
/* Restore the non processed mass (integrated) */ /* Restore the non processed mass (integrated) */
if (snii->integrated.ejected_mass.data != NULL) { if (snii->integrated.ejected_mass_non_processed.data != NULL) {
snii->integrated.ejected_mass.data = snii->integrated.ejected_mass_non_processed.data =
(float *)malloc(sizeof(float) * snii->integrated.ejected_mass.N); (float *)malloc(sizeof(float) * snii->integrated.ejected_mass_non_processed.N);
if (snii->integrated.ejected_mass.data == NULL) { if (snii->integrated.ejected_mass_non_processed.data == NULL) {
error("Failed to allocate memory for the yields"); error("Failed to allocate memory for the yields");
} }
restart_read_blocks((void *)snii->integrated.ejected_mass.data, restart_read_blocks((void *)snii->integrated.ejected_mass_non_processed.data,
sizeof(float), snii->integrated.ejected_mass.N, stream, sizeof(float), snii->integrated.ejected_mass_non_processed.N, stream,
NULL, "non_processed_mass_int"); NULL, "non_processed_mass_int");
} }
/* Restore the non processed mass (raw) */ /* Restore the non processed mass (raw) */
if (snii->raw.ejected_mass.data != NULL) { if (snii->raw.ejected_mass_non_processed.data != NULL) {
snii->raw.ejected_mass.data = snii->raw.ejected_mass_non_processed.data =
(float *)malloc(sizeof(float) * snii->raw.ejected_mass.N); (float *)malloc(sizeof(float) * snii->raw.ejected_mass_non_processed.N);
if (snii->raw.ejected_mass.data == NULL) { if (snii->raw.ejected_mass_non_processed.data == NULL) {
error("Failed to allocate memory for the yields"); error("Failed to allocate memory for the yields");
} }
restart_read_blocks((void *)snii->raw.ejected_mass.data, restart_read_blocks((void *)snii->raw.ejected_mass_non_processed.data,
sizeof(float), snii->raw.ejected_mass.N, stream, sizeof(float), snii->raw.ejected_mass_non_processed.N, stream,
NULL, "non_processed_mass_raw"); NULL, "non_processed_mass_raw");
} }
...@@ -591,6 +591,6 @@ void supernovae_ii_clean(struct supernovae_ii *snii) { ...@@ -591,6 +591,6 @@ void supernovae_ii_clean(struct supernovae_ii *snii) {
interpolate_1d_free(&snii->integrated.ejected_mass_processed); interpolate_1d_free(&snii->integrated.ejected_mass_processed);
interpolate_1d_free(&snii->raw.ejected_mass_processed); interpolate_1d_free(&snii->raw.ejected_mass_processed);
interpolate_1d_free(&snii->integrated.ejected_mass); interpolate_1d_free(&snii->integrated.ejected_mass_non_processed);
interpolate_1d_free(&snii->raw.ejected_mass); interpolate_1d_free(&snii->raw.ejected_mass_non_processed);
} }
...@@ -37,9 +37,9 @@ void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, fl ...@@ -37,9 +37,9 @@ void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, fl
float log_m2, float *yields); float log_m2, float *yields);
void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m, void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m,
float *yields); float *yields);
float supernovae_ii_get_ejected_mass_fraction_from_integral(const struct supernovae_ii *snii, float supernovae_ii_get_ejected_mass_fraction_non_processed_from_integral(const struct supernovae_ii *snii,
float log_m1, float log_m2); float log_m1, float log_m2);
float supernovae_ii_get_ejected_mass_fraction_from_raw(const struct supernovae_ii *snii, float supernovae_ii_get_ejected_mass_fraction_non_processed_from_raw(const struct supernovae_ii *snii,
float log_m); float log_m);
float supernovae_ii_get_ejected_mass_fraction_processed_from_integral( float supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
......
...@@ -59,10 +59,10 @@ INLINE static int star_formation_is_star_forming( ...@@ -59,10 +59,10 @@ INLINE static int star_formation_is_star_forming(
const struct cooling_function_data* restrict cooling, const struct cooling_function_data* restrict cooling,
const struct entropy_floor_properties* restrict entropy_floor) { const struct entropy_floor_properties* restrict entropy_floor) {
/* /\* Check if collapsing particles *\/ */ /* Check if collapsing particles */
/* if (xp->sf_data.div_v > 0) { */ if (xp->sf_data.div_v > 0) {
/* return 0; */ return 0;
/* } */ }
const float temperature = cooling_get_temperature(phys_const, hydro_props, us, const float temperature = cooling_get_temperature(phys_const, hydro_props, us,
cosmo, cooling, p, xp); cosmo, cooling, p, xp);
......
...@@ -47,15 +47,12 @@ INLINE static void stars_read_particles(struct spart *sparts, ...@@ -47,15 +47,12 @@ INLINE static void stars_read_particles(struct spart *sparts,
UNIT_CONV_NO_UNITS, sparts, id); UNIT_CONV_NO_UNITS, sparts, id);
list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL, list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
UNIT_CONV_LENGTH, sparts, h); UNIT_CONV_LENGTH, sparts, h);
// TODO take it from initial mass
list[5] = io_make_input_field("BirthMass", FLOAT, 1, COMPULSORY, list[5] = io_make_input_field("BirthMass", FLOAT, 1, OPTIONAL,
UNIT_CONV_MASS, sparts, sf_data.birth_mass); UNIT_CONV_MASS, sparts, sf_data.birth_mass);
// TODO make it optional
list[6] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS, list[6] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS,
sparts, birth_time); sparts, birth_time);
// TODO read birth density
} }
INLINE static void convert_spart_pos(const struct engine *e, INLINE static void convert_spart_pos(const struct engine *e,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment