feedback.c 9.9 KB
Newer Older
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2018 Loic Hausammann (loic.hausammann@epfl.ch)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

/* Include header */
#include "feedback.h"

/* Local includes */
#include "cosmology.h"
#include "engine.h"
#include "error.h"
#include "feedback_properties.h"
#include "hydro_properties.h"
#include "part.h"
#include "stellar_evolution.h"
#include "units.h"

#include <strings.h>

/**
 * @brief Update the properties of the particle due to a supernovae.
 *
 * @param p The #part to consider.
 * @param xp The #xpart to consider.
Loic Hausammann's avatar
Loic Hausammann committed
40
 * @param e The #engine.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68
 */
void feedback_update_part(struct part* restrict p, struct xpart* restrict xp,
                          const struct engine* restrict e) {

  /* Did the particle receive a supernovae */
  if (xp->feedback_data.delta_mass == 0) return;

  const struct cosmology* cosmo = e->cosmology;

  /* Turn off the cooling */
  xp->cooling_data.time_last_event = e->time;

  /* Update mass */
  const float old_mass = hydro_get_mass(p);
  const float new_mass = old_mass + xp->feedback_data.delta_mass;

  if (xp->feedback_data.delta_mass < 0.) {
    error("Delta mass smaller than 0");
  }

  hydro_set_mass(p, new_mass);

  xp->feedback_data.delta_mass = 0;

  /* Update the density */
  p->rho *= new_mass / old_mass;

  /* Update internal energy */
Loic Hausammann's avatar
Loic Hausammann committed
69 70
  const float u =
      hydro_get_physical_internal_energy(p, xp, cosmo) * old_mass / new_mass;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89
  const float u_new = u + xp->feedback_data.delta_u;

  hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
  hydro_set_drifted_physical_internal_energy(p, cosmo, u_new);

  xp->feedback_data.delta_u = 0.;

  /* Update the velocities */
  for (int i = 0; i < 3; i++) {
    const float dv = xp->feedback_data.delta_p[i] / new_mass;

    xp->v_full[i] += dv;
    p->v[i] += dv;

    xp->feedback_data.delta_p[i] = 0;
  }
}

/**
Loic Hausammann's avatar
Loic Hausammann committed
90
 * @brief Should we do feedback for this star?
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
91 92
 *
 * @param sp The star to consider.
Loic Hausammann's avatar
Loic Hausammann committed
93 94 95 96
 * @param feedback_props The #feedback_props.
 * @param with_cosmology Is the cosmology switch on?
 * @param cosmo The #cosmology.
 * @param time The current time.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118
 */
int feedback_will_do_feedback(const struct spart* sp,
                              const struct feedback_props* feedback_props,
                              const int with_cosmology,
                              const struct cosmology* cosmo,
                              const double time) {

  return (sp->birth_time != -1.);
}

/**
 * @brief Should this particle be doing any feedback-related operation?
 *
 * @param sp The #spart.
 * @param time The current simulation time (Non-cosmological runs).
 * @param cosmo The cosmological model (cosmological runs).
 * @param with_cosmology Are we doing a cosmological run?
 */
int feedback_is_active(const struct spart* sp, const double time,
                       const struct cosmology* cosmo,
                       const int with_cosmology) {

Loic Hausammann's avatar
Loic Hausammann committed
119
  // TODO improve this with estimates for SNII and SNIa
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
  if (sp->birth_time == -1.) return 0;

  if (with_cosmology) {
    return ((double)cosmo->a) > sp->birth_scale_factor;
  } else {
    return time > sp->birth_time;
  }
}

/**
 * @brief Returns the length of time since the particle last did
 * enrichment/feedback.
 *
 * @param sp The #spart.
 * @param with_cosmology Are we running with cosmological time integration on?
 * @param cosmo The cosmological model.
 * @param time The current time (since the Big Bang / start of the run) in
 * internal units.
 * @param dt_star the length of this particle's time-step in internal units.
 * @return The length of the enrichment step in internal units.
 */
double feedback_get_enrichment_timestep(const struct spart* sp,
                                        const int with_cosmology,
                                        const struct cosmology* cosmo,
                                        const double time,
                                        const double dt_star) {
  return dt_star;
}

/**
 * @brief Prepares a s-particle for its feedback interactions
 *
 * @param sp The particle to act upon
 */
void feedback_init_spart(struct spart* sp) {

  sp->feedback_data.enrichment_weight = 0.f;
}

/**
 * @brief Prepares a star's feedback field before computing what
 * needs to be distributed.
 */
void feedback_reset_feedback(struct spart* sp,
                             const struct feedback_props* feedback_props) {

  /* Zero the energy of supernovae */
  sp->feedback_data.energy_ejected = 0;
}

/**
 * @brief Initialises the s-particles feedback props for the first time
 *
 * This function is called only once just after the ICs have been
 * read in to do some conversions.
 *
 * @param sp The particle to act upon.
 * @param feedback_props The properties of the feedback model.
 */
void feedback_first_init_spart(struct spart* sp,
                               const struct feedback_props* feedback_props) {

  feedback_init_spart(sp);

  feedback_reset_feedback(sp, feedback_props);
}

/**
 * @brief Initialises the s-particles feedback props for the first time
 *
 * This function is called only once just after the ICs have been
 * read in to do some conversions.
 *
 * @param sp The particle to act upon.
 * @param feedback_props The properties of the feedback model.
 */
void feedback_prepare_spart(struct spart* sp,
                            const struct feedback_props* feedback_props) {}

/**
 * @brief Evolve the stellar properties of a #spart.
 *
 * This function compute the SN rate and yields before sending
 * this information to a different MPI rank.
 *
 * @param sp The particle to act upon
 * @param feedback_props The #feedback_props structure.
 * @param cosmo The current cosmological model.
 * @param us The unit system.
 * @param phys_const The #phys_const.
 * @param star_age_beg_step The age of the star at the star of the time-step in
 * internal units.
 * @param dt The time-step size of this star in internal units.
 * @param time The physical time in internal units.
 * @param ti_begin The integer time at the beginning of the step.
 * @param with_cosmology Are we running with cosmology on?
 */
void feedback_evolve_spart(struct spart* restrict sp,
                           const struct feedback_props* feedback_props,
                           const struct cosmology* cosmo,
                           const struct unit_system* us,
                           const struct phys_const* phys_const,
                           const double star_age_beg_step, const double dt,
                           const double time, const integertime_t ti_begin,
                           const int with_cosmology) {

#ifdef SWIFT_DEBUG_CHECKS
  if (sp->birth_time == -1.) error("Evolving a star particle that should not!");
Loic Hausammann's avatar
Loic Hausammann committed
228 229 230 231

  if (star_age_beg_step < -1e-6) {
    error("Negative age for a star");
  }
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
232
#endif
Loic Hausammann's avatar
Loic Hausammann committed
233 234
  const double star_age_beg_step_safe =
      star_age_beg_step < 0 ? 0 : star_age_beg_step;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
235 236 237 238 239 240 241 242 243 244

  /* Reset the feedback */
  feedback_reset_feedback(sp, feedback_props);

  /* Add missing h factor */
  const float hi_inv = 1.f / sp->h;
  const float hi_inv_dim = pow_dimension(hi_inv); /* 1/h^d */

  sp->feedback_data.enrichment_weight *= hi_inv_dim;

Loic Hausammann's avatar
Loic Hausammann committed
245 246 247 248 249 250 251 252
  /* Pick the correct table. (if only one table, threshold is < 0) */
  const float metal = chemistry_get_total_metal_mass_fraction_for_feedback(sp);
  const float threshold = feedback_props->metallicity_max_first_stars;

  const struct stellar_model* model =
      metal < threshold ? &feedback_props->stellar_model_first_stars
                        : &feedback_props->stellar_model;

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
253
  /* Compute the stellar evolution */
Loic Hausammann's avatar
Loic Hausammann committed
254 255
  stellar_evolution_evolve_spart(sp, model, cosmo, us, phys_const, ti_begin,
                                 star_age_beg_step_safe, dt);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273

  /* Transform the number of SN to the energy */
  sp->feedback_data.energy_ejected =
      sp->feedback_data.number_sn * feedback_props->energy_per_supernovae;
}

/**
 * @brief Write a feedback struct to the given FILE as a stream of bytes.
 *
 * @param feedback the struct
 * @param stream the file stream
 */
void feedback_struct_dump(const struct feedback_props* feedback, FILE* stream) {

  restart_write_blocks((void*)feedback, sizeof(struct feedback_props), 1,
                       stream, "feedback", "feedback function");

  stellar_evolution_dump(&feedback->stellar_model, stream);
Loic Hausammann's avatar
Loic Hausammann committed
274 275 276
  if (feedback->metallicity_max_first_stars != -1) {
    stellar_evolution_dump(&feedback->stellar_model_first_stars, stream);
  }
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
277 278 279 280 281 282 283 284 285
}

/**
 * @brief Restore a feedback struct from the given FILE as a stream of
 * bytes.
 *
 * @param feedback the struct
 * @param stream the file stream
 */
286
void feedback_struct_restore(struct feedback_props* feedback, FILE* stream) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
287 288 289 290

  restart_read_blocks((void*)feedback, sizeof(struct feedback_props), 1, stream,
                      NULL, "feedback function");

291
  stellar_evolution_restore(&feedback->stellar_model, stream);
Loic Hausammann's avatar
Loic Hausammann committed
292 293 294 295

  if (feedback->metallicity_max_first_stars != -1) {
    stellar_evolution_restore(&feedback->stellar_model_first_stars, stream);
  }
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
296 297 298 299 300 301 302 303 304 305
}

/**
 * @brief Clean the allocated memory.
 *
 * @param feedback the #feedback_props.
 */
void feedback_clean(struct feedback_props* feedback) {

  stellar_evolution_clean(&feedback->stellar_model);
Loic Hausammann's avatar
Loic Hausammann committed
306 307 308
  if (feedback->metallicity_max_first_stars != -1) {
    stellar_evolution_clean(&feedback->stellar_model_first_stars);
  }
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
309
}