diff --git a/src/chemistry/GEAR/chemistry.h b/src/chemistry/GEAR/chemistry.h
index ebeb41b3ba991ced597791e3e0e1c3625aa34b53..b6c5370addeb3e7d7661da8da5096ce8e42ba5f4 100644
--- a/src/chemistry/GEAR/chemistry.h
+++ b/src/chemistry/GEAR/chemistry.h
@@ -21,7 +21,6 @@
 
 /**
  * @file src/chemistry/GEAR/chemistry.h
- * @brief Empty infrastructure for the cases without chemistry function
  */
 
 /* Some standard headers. */
diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c
index 7bd4959d2a6d6c41027f92315438a1ac1daedfa7..34c8498020f4828fc3fc145a9220b61b65e3a054 100644
--- a/src/cooling/grackle/cooling.c
+++ b/src/cooling/grackle/cooling.c
@@ -655,8 +655,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
                       const struct part* restrict p,
                       struct xpart* restrict xp) {
 
-  error("TODO: use energy after adiabatic cooling");
-
   /* set current time */
   code_units units = cooling->units;
 
@@ -711,31 +709,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
   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.
  *
@@ -797,13 +770,6 @@ void cooling_cool_part(const struct phys_const* restrict phys_const,
   /* Calculate the cooling rate */
   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 */
   /* Hydro */
   double u_ad = u_new + hydro_du_dt * dt_therm;
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index f67affac5e7b929c92e77176355ecc9aed584771..55dbb609f10c7da203e5c626bd7db04d1ad29a1c 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -105,9 +105,6 @@ gr_float cooling_time(const struct phys_const* restrict phys_const,
                       const struct cosmology* restrict cosmo,
                       const struct cooling_function_data* restrict cooling,
                       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,
                        const struct unit_system* restrict us,
diff --git a/src/feedback/GEAR/feedback_iact.h b/src/feedback/GEAR/feedback_iact.h
index 93c37cadc2d1c5e132258bfd6506802e1b3a3848..a11c2bf53b45fd601c67c2ee10d72666bd4936bb 100644
--- a/src/feedback/GEAR/feedback_iact.h
+++ b/src/feedback/GEAR/feedback_iact.h
@@ -115,13 +115,11 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
 
   /* Mass received */
   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 dm = m_ej * weight;
   const double new_mass = mj + dm;
 
   /* Energy received */
-  // TODO compute inverse before feedback loop
   const double du = e_sn * weight / new_mass;
 
   xpj->feedback_data.delta_mass += dm;
diff --git a/src/feedback/GEAR/feedback_properties.h b/src/feedback/GEAR/feedback_properties.h
index 4b2d4f207b779779b099c9b61fc7ee0b41e27d1b..df79e6e797e999a72903842e2219704209633775 100644
--- a/src/feedback/GEAR/feedback_properties.h
+++ b/src/feedback/GEAR/feedback_properties.h
@@ -63,8 +63,6 @@ __attribute__((always_inline)) INLINE static void feedback_props_print(
 /**
  * @brief Initialize the global properties of the feedback scheme.
  *
- * By default, takes the values provided by the hydro.
- *
  * @param fp The #feedback_props.
  * @param phys_const The physical constants in the internal unit system.
  * @param us The internal unit system.
diff --git a/src/feedback/GEAR/interpolation.h b/src/feedback/GEAR/interpolation.h
index ed81e838988357f12a10aad692e9978c7b8b3abb..88d6ed30fc6be2377108e7110758a1ab68f32385 100644
--- a/src/feedback/GEAR/interpolation.h
+++ b/src/feedback/GEAR/interpolation.h
@@ -19,6 +19,9 @@
 #ifndef SWIFT_GEAR_INTERPOLATION_H
 #define SWIFT_GEAR_INTERPOLATION_H
 
+/**
+ * @brief Type of boundary condition available.
+ */
 enum interpolate_boundary_condition {
   /* No extrapolation => raise errors */
   boundary_condition_error,
@@ -33,6 +36,9 @@ enum interpolate_boundary_condition {
   boundary_condition_const,
 };
 
+/**
+ * @brief Structure for the interpolation.
+ */
 struct interpolation_1d {
   /* Data to interpolate */
   float *data;
@@ -53,8 +59,6 @@ struct interpolation_1d {
 /**
  * @brief Initialize the #interpolation_1d.
  *
- * Assumes x are linear in log.
- *
  * @params interp The #interpolation_1d.
  * @params xmin Minimal value of x (in log).
  * @params xmax Maximal value of x (in log).
diff --git a/src/feedback/GEAR/stellar_evolution.c b/src/feedback/GEAR/stellar_evolution.c
index c42baae7ed315ab6a586167f54311ef5f322d447..757bbef2d3a7bcae68fcd76de99312b467ff149d 100644
--- a/src/feedback/GEAR/stellar_evolution.c
+++ b/src/feedback/GEAR/stellar_evolution.c
@@ -111,7 +111,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
 
   /* SNII */
   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 */
   sp->feedback_data.mass_ejected = mass_frac_snii * sp->sf_data.birth_mass
@@ -136,7 +136,7 @@ void stellar_evolution_compute_continuous_feedback_properties(
                            snii_yields);
 
   /* 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);
 
   /* Set the yields */
@@ -222,7 +222,7 @@ void stellar_evolution_compute_discrete_feedback_properties(
 
   /* Compute the mass fraction of non processed elements */
   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);
 
   /* Set the yields */
diff --git a/src/feedback/GEAR/stellar_evolution_struct.h b/src/feedback/GEAR/stellar_evolution_struct.h
index c80e57711f5aef9e65a8726c9911d3ada6096889..3a88b5669d24860ac8361d0a567ee3c823ccc265 100644
--- a/src/feedback/GEAR/stellar_evolution_struct.h
+++ b/src/feedback/GEAR/stellar_evolution_struct.h
@@ -117,12 +117,12 @@ struct supernovae_ii {
   struct {
     /*! Mass fraction of metals ejected by a supernovae. */
     struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-    
+
     /*! Total mass fraction ejected. */
     struct interpolation_1d ejected_mass_processed;
 
     /*! Mass fraction ejected and not processed (=> with the star metallicity). */
-    struct interpolation_1d ejected_mass;
+    struct interpolation_1d ejected_mass_non_processed;
   } raw;
 
   /*! Yields integrated */
@@ -130,12 +130,12 @@ struct supernovae_ii {
     /*! Integrated (over the IMF) mass fraction of metals ejected by a supernovae
      */
     struct interpolation_1d yields[GEAR_CHEMISTRY_ELEMENT_COUNT];
-    
+
     /*! Total mass fraction ejected (integrated over the IMF) */
     struct interpolation_1d ejected_mass_processed;
 
     /*! Mass fraction ejected and not processed (=> with the star metallicity) */
-    struct interpolation_1d ejected_mass;
+    struct interpolation_1d ejected_mass_non_processed;
   } integrated;
 
   /*! Minimal mass for a SNII */
diff --git a/src/feedback/GEAR/supernovae_ii.c b/src/feedback/GEAR/supernovae_ii.c
index 4790a590b41a095dbbdf0facead0777f7bcd8de1..d34ac6dc9f2a1ecb73138597f7b82ae6125739b7 100644
--- a/src/feedback/GEAR/supernovae_ii.c
+++ b/src/feedback/GEAR/supernovae_ii.c
@@ -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 m1 The lower 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 mass_ejected_1 = interpolate_1d(&snii->integrated.ejected_mass, log_m1);
-  float mass_ejected_2 = interpolate_1d(&snii->integrated.ejected_mass, log_m2);
+  float mass_ejected_1 = interpolate_1d(&snii->integrated.ejected_mass_non_processed, log_m1);
+  float mass_ejected_2 = interpolate_1d(&snii->integrated.ejected_mass_non_processed, log_m2);
 
   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 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) {
 
-  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,
       "Ej", &previous_count, interpolation_size);
 
   /* Read the mass ejected of non processed gas */
-  supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass, 
-				  &snii->integrated.ejected_mass,
+  supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass_non_processed,
+				  &snii->integrated.ejected_mass_non_processed,
 				  phys_const, sm, group_id, "Ejnp",
                                   &previous_count, interpolation_size);
 
@@ -449,16 +449,16 @@ void supernovae_ii_dump(const struct supernovae_ii *snii, FILE *stream,
   }
 
   /*! Dump the non processed mass (integrated). */
-  if (snii->integrated.ejected_mass.data != NULL) {
-    restart_write_blocks((void *)snii->integrated.ejected_mass.data,
-                         sizeof(float), snii->integrated.ejected_mass.N, stream,
+  if (snii->integrated.ejected_mass_non_processed.data != NULL) {
+    restart_write_blocks((void *)snii->integrated.ejected_mass_non_processed.data,
+                         sizeof(float), snii->integrated.ejected_mass_non_processed.N, stream,
                          "non_processed_mass_int", "non_processed_mass_int");
   }
 
   /*! Dump the non processed mass (raw). */
-  if (snii->raw.ejected_mass.data != NULL) {
-    restart_write_blocks((void *)snii->raw.ejected_mass.data,
-                         sizeof(float), snii->raw.ejected_mass.N, stream,
+  if (snii->raw.ejected_mass_non_processed.data != NULL) {
+    restart_write_blocks((void *)snii->raw.ejected_mass_non_processed.data,
+                         sizeof(float), snii->raw.ejected_mass_non_processed.N, stream,
                          "non_processed_mass_raw", "non_processed_mass_raw");
   }
 
@@ -550,28 +550,28 @@ void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream,
   }
 
   /* Restore the non processed mass (integrated) */
-  if (snii->integrated.ejected_mass.data != NULL) {
-    snii->integrated.ejected_mass.data =
-        (float *)malloc(sizeof(float) * snii->integrated.ejected_mass.N);
-      if (snii->integrated.ejected_mass.data == NULL) {
+  if (snii->integrated.ejected_mass_non_processed.data != NULL) {
+    snii->integrated.ejected_mass_non_processed.data =
+      (float *)malloc(sizeof(float) * snii->integrated.ejected_mass_non_processed.N);
+    if (snii->integrated.ejected_mass_non_processed.data == NULL) {
 	error("Failed to allocate memory for the yields");
       }
 
-    restart_read_blocks((void *)snii->integrated.ejected_mass.data,
-                        sizeof(float), snii->integrated.ejected_mass.N, stream,
+    restart_read_blocks((void *)snii->integrated.ejected_mass_non_processed.data,
+                        sizeof(float), snii->integrated.ejected_mass_non_processed.N, stream,
                         NULL, "non_processed_mass_int");
   }
 
   /* Restore the non processed mass (raw) */
-  if (snii->raw.ejected_mass.data != NULL) {
-    snii->raw.ejected_mass.data =
-        (float *)malloc(sizeof(float) * snii->raw.ejected_mass.N);
-      if (snii->raw.ejected_mass.data == NULL) {
+  if (snii->raw.ejected_mass_non_processed.data != NULL) {
+    snii->raw.ejected_mass_non_processed.data =
+      (float *)malloc(sizeof(float) * snii->raw.ejected_mass_non_processed.N);
+    if (snii->raw.ejected_mass_non_processed.data == NULL) {
 	error("Failed to allocate memory for the yields");
       }
 
-    restart_read_blocks((void *)snii->raw.ejected_mass.data,
-                        sizeof(float), snii->raw.ejected_mass.N, stream,
+    restart_read_blocks((void *)snii->raw.ejected_mass_non_processed.data,
+                        sizeof(float), snii->raw.ejected_mass_non_processed.N, stream,
                         NULL, "non_processed_mass_raw");
   }
 
@@ -591,6 +591,6 @@ void supernovae_ii_clean(struct supernovae_ii *snii) {
 
   interpolate_1d_free(&snii->integrated.ejected_mass_processed);
   interpolate_1d_free(&snii->raw.ejected_mass_processed);
-  interpolate_1d_free(&snii->integrated.ejected_mass);
-  interpolate_1d_free(&snii->raw.ejected_mass);
+  interpolate_1d_free(&snii->integrated.ejected_mass_non_processed);
+  interpolate_1d_free(&snii->raw.ejected_mass_non_processed);
 }
diff --git a/src/feedback/GEAR/supernovae_ii.h b/src/feedback/GEAR/supernovae_ii.h
index 3f525df502e44a6185ce9837ec665fb716651982..278433010365c27ea9130b82e47e7849968aaf5c 100644
--- a/src/feedback/GEAR/supernovae_ii.h
+++ b/src/feedback/GEAR/supernovae_ii.h
@@ -37,9 +37,9 @@ void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii, fl
 					    float log_m2, float *yields);
 void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii, float log_m,
 				       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 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 supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
diff --git a/src/star_formation/GEAR/star_formation.h b/src/star_formation/GEAR/star_formation.h
index 10ab45e50db84e216d9226704640fe39a73c32ea..82632dcb06043a42de4f4cf6a8b84a078492997a 100644
--- a/src/star_formation/GEAR/star_formation.h
+++ b/src/star_formation/GEAR/star_formation.h
@@ -59,10 +59,10 @@ INLINE static int star_formation_is_star_forming(
     const struct cooling_function_data* restrict cooling,
     const struct entropy_floor_properties* restrict entropy_floor) {
 
-  /* /\* Check if collapsing particles *\/ */
-  /* if (xp->sf_data.div_v > 0) { */
-  /*   return 0; */
-  /* } */
+  /* Check if collapsing particles */
+  if (xp->sf_data.div_v > 0) {
+    return 0;
+  }
 
   const float temperature = cooling_get_temperature(phys_const, hydro_props, us,
                                                     cosmo, cooling, p, xp);
diff --git a/src/stars/GEAR/stars_io.h b/src/stars/GEAR/stars_io.h
index d5ad49f953b2f49b92a54c63168689328b06d224..27a7a8f241ecc8536d1c9ed8016e4bed29762143 100644
--- a/src/stars/GEAR/stars_io.h
+++ b/src/stars/GEAR/stars_io.h
@@ -47,15 +47,12 @@ INLINE static void stars_read_particles(struct spart *sparts,
                                 UNIT_CONV_NO_UNITS, sparts, id);
   list[4] = io_make_input_field("SmoothingLength", FLOAT, 1, OPTIONAL,
                                 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);
 
-  // TODO make it optional
   list[6] = io_make_input_field("BirthTime", FLOAT, 1, OPTIONAL, UNIT_CONV_MASS,
                                 sparts, birth_time);
-
-  // TODO read birth density
 }
 
 INLINE static void convert_spart_pos(const struct engine *e,