GEAR: fixing physics

src/cooling/grackle/cooling.c

@@ -167,7 +167,7 @@ void cooling_compute_equilibrium(
   double dt = fabs(cooling_time(phys_const, us, hydro_properties, cosmo,
                                 &cooling_tmp, &p_tmp, xp));
   cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp,
-                     &p_tmp, xp, dt);
+                     &p_tmp, xp, dt, dt);
   dt = alpha * fabs(cooling_time(phys_const, us, hydro_properties, cosmo,
                                  &cooling_tmp, &p_tmp, xp));
 
@@ -184,7 +184,7 @@ void cooling_compute_equilibrium(
 
     /* update chemistry */
     cooling_new_energy(phys_const, us, cosmo, hydro_properties, &cooling_tmp,
-                       &p_tmp, xp, dt);
+                       &p_tmp, xp, dt, dt);
   } while (step < max_step && !cooling_converged(xp, &old, conv_limit));
 
   if (step == max_step)
@@ -572,6 +572,7 @@ void cooling_apply_self_shielding(
  * @param p Pointer to the particle data.
  * @param xp Pointer to the particle extra data
  * @param dt The time-step of this particle.
+ * @param dt_therm The time-step operator used for thermal quantities.
  *
  * @return du / dt
  */
@@ -581,7 +582,8 @@ gr_float cooling_new_energy(
     const struct cosmology* restrict cosmo,
     const struct hydro_props* hydro_props,
     const struct cooling_function_data* restrict cooling,
-    const struct part* restrict p, struct xpart* restrict xp, double dt) {
+    const struct part* restrict p, struct xpart* restrict xp, double dt,
+    double dt_therm) {
 
   /* set current time */
   code_units units = cooling->units;
@@ -605,11 +607,7 @@ gr_float cooling_new_energy(
   /* general particle data */
   gr_float density = hydro_get_physical_density(p, cosmo);
   gr_float energy = hydro_get_physical_internal_energy(p, xp, cosmo) +
-                    dt * hydro_get_physical_internal_energy_dt(p, cosmo);
-
-  /* We now need to check that we are not going to go below any of the limits */
-  const double u_minimal = hydro_props->minimal_internal_energy;
-  energy = max(energy, u_minimal);
+                    dt_therm * hydro_get_physical_internal_energy_dt(p, cosmo);
 
   /* initialize density */
   data.density = &density;
@@ -771,33 +769,63 @@ void cooling_cool_part(const struct phys_const* restrict phys_const,
     return;
   }
 
-  /* Get the change in internal energy due to hydro forces */
-  const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
-
   /* Current energy */
   const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo);
 
+  /* Energy after the adiabatic cooling */
+  float u_ad_before = u_old +
+    dt_therm * hydro_get_physical_internal_energy_dt(p, cosmo);
+
+  /* We now need to check that we are not going to go below any of the limits */
+  const double u_minimal = hydro_props->minimal_internal_energy;
+  if (u_ad_before < u_minimal) {
+    u_ad_before = u_minimal;
+    const float du_dt = (u_ad_before - u_old) / dt_therm;
+    hydro_set_physical_internal_energy_dt(p, cosmo, du_dt);
+  }
+
   /* Calculate energy after dt */
   gr_float u_new = cooling_new_energy(phys_const, us, cosmo, hydro_props,
-                                      cooling, p, xp, dt);
+                                      cooling, p, xp, dt, dt_therm);
+
+  /* Get the change in internal energy due to hydro forces */
+  float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
 
   /* We now need to check that we are not going to go below any of the limits */
-  const double u_minimal = hydro_props->minimal_internal_energy;
   u_new = max(u_new, u_minimal);
 
-  /* Expected change in energy over the next kick step
-     (assuming no change in dt) */
-  const double delta_u = u_new - u_old;
+  /* 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;
+  if (u_ad < u_minimal) {
+    hydro_du_dt = (u_ad - u_new) / dt_therm;
+    u_ad = u_minimal;
+  }
+
+  /* Cooling */
+  const double u_cool = u_ad + cool_du_dt * dt_therm;
+  if (u_cool < u_minimal) {
+    cool_du_dt = (u_cool - u_ad) / dt_therm;
+  }
 
-  /* Turn this into a rate of change (including cosmology term) */
-  const float cooling_du_dt = delta_u / dt_therm;
+  cool_du_dt = (u_new - u_old) / dt_therm;
 
   /* Update the internal energy time derivative */
-  cooling_apply(p, xp, cosmo, cooling_du_dt, u_new);
+  hydro_set_physical_internal_energy_dt(p, cosmo, cool_du_dt /* + hydro_du_dt */);
 
   /* Store the radiated energy */
   xp->cooling_data.radiated_energy -=
-      hydro_get_mass(p) * (cooling_du_dt - hydro_du_dt) * dt;
+      hydro_get_mass(p) * cool_du_dt * dt_therm;
 }
 
 /**

src/cooling/grackle/cooling.h

@@ -96,7 +96,8 @@ gr_float cooling_new_energy(
     const struct cosmology* restrict cosmo,
     const struct hydro_props* hydro_properties,
     const struct cooling_function_data* restrict cooling,
-    const struct part* restrict p, struct xpart* restrict xp, double dt);
+    const struct part* restrict p, struct xpart* restrict xp, double dt,
+    double dt_therm);
 
 gr_float cooling_time(const struct phys_const* restrict phys_const,
                       const struct unit_system* restrict us,

src/feedback/GEAR/feedback.c

@@ -111,6 +111,7 @@ int feedback_is_active(const struct spart* sp, const double time,
                        const struct cosmology* cosmo,
                        const int with_cosmology) {
 
+  // TODO improve this with estimates for SNII and SNIa
   if (sp->birth_time == -1.) return 0;
 
   if (with_cosmology) {

src/feedback/GEAR/initial_mass_function.c

@@ -75,19 +75,22 @@ void initial_mass_function_print(const struct initial_mass_function *imf) {
  * The x are supposed to be linear in log.
  *
  * @param imf The #initial_mass_function.
- * @param interp The #interpolation_1d.
+ * @param data The data to integrate.
+ * @param count The number of element in data.
+ * @param log_mass_min The value of the first element.
+ * @param step_size The distance between two points.
  */
 void initial_mass_function_integrate(const struct initial_mass_function *imf,
-                                     struct interpolation_1d *interp) {
+                                     float *data, size_t count, float log_mass_min, float step_size) {
 
   /* Index in the data */
-  int j = 1;
-  const float mass_min = pow(10, interp->xmin);
-  const float mass_max = pow(10, interp->xmin + (interp->N - 1) * interp->dx);
+  size_t j = 1;
+  const float mass_min = pow(10, log_mass_min);
+  const float mass_max = pow(10, log_mass_min + (count - 1) * step_size);
 
   float m = mass_min;
 
-  float *tmp = (float *)malloc(sizeof(float) * interp->N);
+  float *tmp = (float *)malloc(sizeof(float) * count);
 
   /* Set lower limit */
   tmp[0] = 0;
@@ -104,12 +107,12 @@ void initial_mass_function_integrate(const struct initial_mass_function *imf,
     }
 
     /* Integrate the data */
-    while (m < imf->mass_limits[i + 1] && j < interp->N) {
+    while (m < imf->mass_limits[i + 1] && j < count) {
 
       /* Compute the masses */
-      const float log_m1 = interp->xmin + (j - 1) * interp->dx;
+      const float log_m1 = log_mass_min + (j - 1) * step_size;
       const float m1 = pow(10, log_m1);
-      const float log_m2 = interp->xmin + j * interp->dx;
+      const float log_m2 = log_mass_min + j * step_size;
       const float m2 = pow(10, log_m2);
       const float dm = m2 - m1;
       const float imf_1 = imf->coef[i] * pow(m1, imf->exp[i]);
@@ -125,7 +128,7 @@ void initial_mass_function_integrate(const struct initial_mass_function *imf,
       /* Compute the integral */
       tmp[j] =
           tmp[j - 1] +
-          0.5 * (imf_1 * interp->data[j - 1] + imf_2 * interp->data[j]) * dm;
+          0.5 * (imf_1 * data[j - 1] + imf_2 * data[j]) * dm;
 
       /* Update j and m */
       j += 1;
@@ -134,15 +137,12 @@ void initial_mass_function_integrate(const struct initial_mass_function *imf,
   }
 
   /* The rest is extrapolated with 0 */
-  for (int k = j; k < interp->N; k++) {
+  for (size_t k = j; k < count; k++) {
     tmp[k] = tmp[k - 1];
   }
 
   /* Copy temporary array */
-  memcpy(interp->data, tmp, interp->N * sizeof(float));
-
-  /* Update the boundary conditions */
-  interp->boundary_condition = boundary_condition_zero_const;
+  memcpy(data, tmp, count * sizeof(float));
 
   /* clean everything */
   free(tmp);

src/feedback/GEAR/initial_mass_function.h

@@ -27,7 +27,7 @@ float initial_mass_function_get_exponent(
 void initial_mass_function_print(const struct initial_mass_function *imf);
 
 void initial_mass_function_integrate(const struct initial_mass_function *imf,
-                                     struct interpolation_1d *interp);
+                                     float *data, size_t count, float log_mass_min, float step_size);
 float initial_mass_function_get_coefficient(
     const struct initial_mass_function *imf, float mass_min, float mass_max);
 float initial_mass_function_get_integral_xi(

src/feedback/GEAR/stellar_evolution.c

@@ -106,19 +106,15 @@ void stellar_evolution_compute_continuous_feedback_properties(
 
   /* Compute the mass ejected */
   /* SNIa */
-  const float mass_frac_snia =
-      supernovae_ia_get_ejected_mass_processed(&sm->snia) *
-      supernovae_ia_get_number(&sm->snia, m_end_step, m_beg_step);
+  const float mass_snia =
+    supernovae_ia_get_ejected_mass_processed(&sm->snia) * number_snia;
 
   /* SNII */
-  const float mass_frac_snii =
-      supernovae_ii_get_ejected_mass_fraction_processed(
-          &sm->snii, log_m_end_step, log_m_beg_step);
-
-  sp->feedback_data.mass_ejected = (mass_frac_snia + mass_frac_snii) * m_init;
+  const float mass_frac_snii = supernovae_ii_get_ejected_mass_fraction_processed(
+          &sm->snii, log_m_end_step, log_m_beg_step); 
 
-  /* Transform into internal units */
-  sp->feedback_data.mass_ejected *= phys_const->const_solar_mass;
+  /* Sum the contributions from SNIa and SNII */
+  sp->feedback_data.mass_ejected = mass_frac_snii * sp->birth.mass + mass_snia;
 
   if (sp->mass <= sp->feedback_data.mass_ejected) {
     error("Stars cannot have negative mass. (%g <= %g). Initial mass = %g",
@@ -146,15 +142,20 @@ void stellar_evolution_compute_continuous_feedback_properties(
   for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
     /* Compute the mass fraction of metals */
     sp->feedback_data.metal_mass_ejected[i] =
-        /* Supernovae Ia yields */
-        snia_yields[i] * number_snia +
-        /* Supernovae II yields */
-        snii_yields[i] +
-        /* Gas contained in stars initial metallicity */
-        chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
+      /* Supernovae II yields */
+      snii_yields[i] +
+      /* Gas contained in stars initial metallicity */
+      chemistry_get_metal_mass_fraction_for_feedback(sp)[i] * non_processed;
 
     /* Convert it to total mass */
     sp->feedback_data.metal_mass_ejected[i] *= sp->sf_data.birth_mass;
+
+    sp->feedback_data.metal_mass_ejected[i] *= sp->sf_data.birth_mass;
+
+    /* Add the Supernovae Ia */
+    sp->feedback_data.metal_mass_ejected[i] +=
+      snia_yields[i] * number_snia * phys_const->const_solar_mass;
+
   }
 }
 
@@ -187,7 +188,7 @@ void stellar_evolution_compute_discrete_feedback_properties(
       number_snii == 0
           ? 0.
           : number_snii /
-                (supernovae_ii_get_number(&sm->snii, m_end_step, m_beg_step) *
+                (supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) *
                  m_init);
 
   /* Compute the mass ejected */
@@ -308,11 +309,15 @@ void stellar_evolution_evolve_spart(
   if (can_produce_snia) {
     /* Compute rates */
     const float number_snia_f =
-        supernovae_ia_get_number(&sm->snia, m_end_step, m_beg_step) * m_init;
+        supernovae_ia_get_number_per_unit_mass(&sm->snia, m_end_step, m_beg_step) * m_init;
 
     /* Get the integer number of supernovae */
     number_snia = stellar_evolution_compute_integer_number_supernovae(
         sp, number_snia_f, ti_begin, random_number_stellar_feedback_1);
+
+    if (number_snia != 0) {
+      message("Exploding %i SNIa", number_snia);
+    }
   }
 
   /* Compute number of SNII */
@@ -320,11 +325,14 @@ void stellar_evolution_evolve_spart(
   if (can_produce_snii) {
     /* Compute rates */
     const float number_snii_f =
-        supernovae_ii_get_number(&sm->snii, m_end_step, m_beg_step) * m_init;
+        supernovae_ii_get_number_per_unit_mass(&sm->snii, m_end_step, m_beg_step) * m_init;
 
     /* Get the integer number of supernovae */
     number_snii = stellar_evolution_compute_integer_number_supernovae(
         sp, number_snii_f, ti_begin, random_number_stellar_feedback_2);
+    if (number_snii != 0) {
+      message("Exploding %i SNII", number_snii);
+    }
   }
 
   /* Does this star produce a supernovae? */

src/feedback/GEAR/supernovae_ia.c

@@ -110,7 +110,7 @@ float supernovae_ia_get_companion_fraction(const struct supernovae_ia *snia,
 
   const float tmp =
       pow(m2, snia->companion_exponent) - pow(m1, snia->companion_exponent);
-  return snia->companion[companion_type].coef * tmp / snia->companion_exponent;
+  return snia->companion[companion_type].coef * tmp;
 }
 
 /**
@@ -123,8 +123,8 @@ float supernovae_ia_get_companion_fraction(const struct supernovae_ia *snia,
  *
  * @return The number of supernovae Ia per unit of mass.
  */
-float supernovae_ia_get_number(const struct supernovae_ia *snia, float m1,
-                               float m2) {
+float supernovae_ia_get_number_per_unit_mass(const struct supernovae_ia *snia, float m1,
+					     float m2) {
 
 #ifdef SWIFT_DEBUG_CHECKS
   if (m1 > m2) error("Mass 1 larger than mass 2 %g > %g.", m1, m2);
@@ -220,13 +220,18 @@ void supernovae_ia_read_yields(struct supernovae_ia *snia,
  */
 void supernovae_ia_init_companion(struct supernovae_ia *snia) {
 
+  /* Get the exponent for the integral */
+  const float exp = snia->companion_exponent;
+  
   for (int i = 0; i < GEAR_NUMBER_TYPE_OF_COMPANION; i++) {
+
     /* Compute the integral */
-    float integral = supernovae_ia_get_companion_fraction(
-        snia, snia->companion[i].mass_min, snia->companion[i].mass_max, i);
+    float integral = pow(snia->companion[i].mass_max, exp + 1.);
+    integral -= pow(snia->companion[i].mass_min, exp + 1.);
+    integral /= exp + 1.;
 
     /* Update the coefficient for a normalization to 1 of the IMF */
-    snia->companion[i].coef *= snia->companion[i].coef / integral;
+    snia->companion[i].coef /=  exp * integral;
   }
 }
 

src/feedback/GEAR/supernovae_ia.h

@@ -32,8 +32,8 @@ float supernovae_ia_get_ejected_mass_processed(
 float supernovae_ia_get_companion_fraction(const struct supernovae_ia *snia,
                                            float m1, float m2,
                                            int companion_type);
-float supernovae_ia_get_number(const struct supernovae_ia *snia, float m1,
-                               float m2);
+float supernovae_ia_get_number_per_unit_mass(const struct supernovae_ia *snia, float m1,
+					     float m2);
 
 void supernovae_ia_read_yields(struct supernovae_ia *snia,
                                struct swift_params *params,

src/feedback/GEAR/supernovae_ii.c

@@ -68,8 +68,8 @@ int supernovae_ii_can_explode(const struct supernovae_ii *snii, float m_low,
  *
  * @return The number of supernovae II per unit of mass.
  */
-float supernovae_ii_get_number(const struct supernovae_ii *snii, float m1,
-                               float m2) {
+float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii, float m1,
+					     float m2) {
 #ifdef SWIFT_DEBUG_CHECKS
   if (m1 > m2) error("Mass 1 larger than mass 2 %g > %g.", m1, m2);
 #endif
@@ -191,13 +191,15 @@ void supernovae_ii_read_yields_array(
   /* Read the dataset */
   io_read_array_dataset(group_id, hdf5_dataset_name, FLOAT, data, count);
 
+  /* Integrate the yields */
+  initial_mass_function_integrate(&sm->imf, data, count, log_mass_min, step_size);
+  // TODO: decrease count in order to keep the same distance between points
+
   /* Initialize the interpolation */
   interpolate_1d_init(interp, log10(snii->mass_min), log10(snii->mass_max),
                       interpolation_size, log_mass_min, step_size, count, data,
-                      boundary_condition_zero);
+		      boundary_condition_zero_const);
 
-  /* Integrate the yields */
-  initial_mass_function_integrate(&sm->imf, interp);
 
   /* Cleanup the memory */
   free(data);

src/feedback/GEAR/supernovae_ii.h

@@ -31,8 +31,8 @@
 void supernovae_ii_print(const struct supernovae_ii *snii);
 int supernovae_ii_can_explode(const struct supernovae_ii *snii, float m_low,
                               float m_high);
-float supernovae_ii_get_number(const struct supernovae_ii *snii, float m1,
-                               float m2);
+float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii, float m1,
+					     float m2);
 void supernovae_ii_get_yields(const struct supernovae_ii *snii, float log_m1,
                               float log_m2, float *yields);
 float supernovae_ii_get_ejected_mass_fraction(const struct supernovae_ii *snii,

src/pressure_floor/GEAR/pressure_floor.h

@@ -75,8 +75,8 @@ pressure_floor_get_physical_pressure(const struct part* p,
   const float rho = hydro_get_physical_density(p, cosmo);
 
   /* Compute the pressure floor */
-  float floor = H_phys * H_phys * rho * pressure_floor_props.constants -
-                p->pressure_floor_data.sigma2;
+  float floor = H_phys * H_phys * rho * pressure_floor_props.constants;// -
+  //p->pressure_floor_data.sigma2;
   floor *= rho * hydro_one_over_gamma;
 
   return fmaxf(pressure_physical, floor);
@@ -98,13 +98,13 @@ pressure_floor_get_comoving_pressure(const struct part* p,
                                      const float pressure_comoving,
                                      const struct cosmology* cosmo) {
 
-  const float a_coef = pow_three_gamma_minus_one(cosmo->a);
+  const float a_coef = pow_three_gamma_minus_one(cosmo->a) * cosmo->a_inv;
   const float rho = hydro_get_comoving_density(p);
 
   /* Compute the pressure floor */
   float floor = kernel_gamma * kernel_gamma * p->h * p->h * rho *
                 pressure_floor_props.constants;
-  floor -= p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
+  //floor -= p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
   floor *= a_coef * rho * hydro_one_over_gamma;
 
   return fmaxf(pressure_comoving, floor);

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);
@@ -75,10 +75,10 @@ INLINE static int star_formation_is_star_forming(
   }
 
   /* Get the required variables */
-  const float sigma2 = p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
+  //const float sigma2 = p->pressure_floor_data.sigma2 * cosmo->a * cosmo->a;
   const float n_jeans_2_3 = starform->n_jeans_2_3;
 
-  const float h = p->h * kernel_gamma;
+  const float h = p->h * kernel_gamma * cosmo->a;
   const float density = hydro_get_physical_density(p, cosmo);
 
   // TODO use GRACKLE */
@@ -89,8 +89,8 @@ INLINE static int star_formation_is_star_forming(
       M_PI_4 / (phys_const->const_newton_G * n_jeans_2_3 * h * h);
   const float density_criterion =
       coef * (hydro_gamma * phys_const->const_boltzmann_k * temperature /
-                  (mu * phys_const->const_proton_mass) +
-              sigma2);
+	      (mu * phys_const->const_proton_mass));// +
+  //sigma2);
 
   /* Check the density criterion */
   return density > density_criterion;