supernovae_ia.c 11.6 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 40 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 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2019 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 "supernovae_ia.h"

/* Local headers */
#include "hdf5_functions.h"
#include "stellar_evolution.h"
#include "stellar_evolution_struct.h"

/**
 * @brief Print the supernovae Ia model.
 *
 * @param snia The #supernovae_ia.
 */
void supernovae_ia_print(const struct supernovae_ia *snia) {

  /* Only the master print */
  if (engine_rank != 0) {
    return;
  }

  message("Mass of the white dwarf = %g", snia->mass_white_dwarf);
  message("Mass range of the progenitor = [%g, %g]", snia->mass_min_progenitor,
          snia->mass_max_progenitor);

  for (int i = 0; i < GEAR_NUMBER_TYPE_OF_COMPANION; i++) {
    message("Mass range of the companion %i = [%g, %g]", i,
            snia->companion[i].mass_min, snia->companion[i].mass_max);
  }
}

/**
 * @brief Check if the given mass is able to produce a SNIa.
 *
 * @param snia The #supernovae_ia model.
 * @param m_low The lower mass.
 * @param m_high The higher mass.
 *
 * @return If the mass is in the range of SNIa.
 */
int supernovae_ia_can_explode(const struct supernovae_ia *snia, float m_low,
                              float m_high) {

  if (m_low > snia->mass_max_progenitor) return 0;

  for (int i = 0; i < GEAR_NUMBER_TYPE_OF_COMPANION; i++) {
    if (m_low < snia->companion[i].mass_max &&
        m_high > snia->companion[i].mass_min) {
      return 1;
    }
  }

  return 0;
}

/**
 * @brief Get the yields of a supernovae Ia.
 *
 * @param snia The #supernovae_ia model.
 */
const float *supernovae_ia_get_yields(const struct supernovae_ia *snia) {
  return snia->yields;
}

/**
 * @brief Get the processed mass ejected of a supernovae Ia.
 *
 * @param snia The #supernovae_ia model.
 */
float supernovae_ia_get_ejected_mass_processed(
    const struct supernovae_ia *snia) {
  return snia->mass_white_dwarf;
}

/**
 * @brief Compute the companion integral (second integral in equation 3.46 in
 * Poirier 2004)
 *
 * @param snia The #supernovae_ia model.
 * @param m1 The lower mass limit.
 * @param m2 The upper mass limit.
 * @param companion_type The type of companion (e.g. index of snia->companion).
 *
 * @return The fraction of companion.
 */
float supernovae_ia_get_companion_fraction(const struct supernovae_ia *snia,
                                           float m1, float m2,
                                           int companion_type) {
#ifdef SWIFT_DEBUG_CHECKS
  if (m1 > m2) error("Mass 1 larger than mass 2 %g > %g.", m1, m2);
#endif

  const float tmp =
      pow(m2, snia->companion_exponent) - pow(m1, snia->companion_exponent);
Loic Hausammann's avatar
Loic Hausammann committed
113
  return snia->companion[companion_type].coef * tmp;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
114 115 116 117 118 119 120 121 122 123 124 125
}

/**
 * @brief Compute the number of supernovae Ia per unit of mass (equation 3.46 in
 * Poirier 2004).
 *
 * @param snia The #supernovae_ia model.
 * @param m1 The lower mass limit.
 * @param m2 The upper mass limit.
 *
 * @return The number of supernovae Ia per unit of mass.
 */
Loic Hausammann's avatar
Loic Hausammann committed
126 127
float supernovae_ia_get_number_per_unit_mass(const struct supernovae_ia *snia,
                                             float m1, float m2) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
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

#ifdef SWIFT_DEBUG_CHECKS
  if (m1 > m2) error("Mass 1 larger than mass 2 %g > %g.", m1, m2);
#endif

  /* Do we have white dwarf? */
  if (m1 > snia->mass_max_progenitor) {
    return 0.;
  }

  float number_companion = 0.;
  for (int i = 0; i < GEAR_NUMBER_TYPE_OF_COMPANION; i++) {
    /* Check if we are in the possible interval */
    if (m1 > snia->companion[i].mass_max || m2 < snia->companion[i].mass_min)
      continue;

    /* Get mass limits */
    const float mass_min = max(m1, snia->companion[i].mass_min);
    const float mass_max = min(m2, snia->companion[i].mass_max);

    /* Compute number of companions */
    number_companion +=
        supernovae_ia_get_companion_fraction(snia, mass_min, mass_max, i);
  }

  /* Use only the white dwarf already created */
  const float mass_min = max(m1, snia->mass_min_progenitor);

  /* Compute number of white dwarf */
  float number_white_dwarf =
      pow(snia->mass_max_progenitor, snia->progenitor_exponent);
  number_white_dwarf -= pow(mass_min, snia->progenitor_exponent);
  number_white_dwarf *= snia->progenitor_coef_exp;

  return number_companion * number_white_dwarf;
};

/**
 * @brief Read the SNIa yields from the table.
 *
 * @param snia The #supernovae_ia model.
 * @param params The #swift_params.
 * @param sm The #stellar_model.
Loic Hausammann's avatar
Loic Hausammann committed
171
 * @param filename The filename of the chemistry table.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
172 173 174
 */
void supernovae_ia_read_yields(struct supernovae_ia *snia,
                               struct swift_params *params,
Loic Hausammann's avatar
Loic Hausammann committed
175 176
                               const struct stellar_model *sm,
                               const char *filename) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
177 178 179 180

  hid_t file_id, group_id;

  /* Open IMF group */
181
  h5_open_group(filename, "Data/SNIa/Metals", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
182

Loic Hausammann's avatar
Loic Hausammann committed
183 184 185 186 187 188 189
  /* Get the number of elements. */
  const hid_t attr = H5Aopen(group_id, "elts", H5P_DEFAULT);
  if (attr < 0) error("Error while opening attribute elts in Data/SNIa/Metals");

  size_t nval = io_get_number_element_in_attribute(attr);
  H5Aclose(attr);

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
190
  /* Read the yields */
Loic Hausammann's avatar
Loic Hausammann committed
191 192
  float *yields = (float *)malloc(sizeof(float) * nval);
  io_read_array_attribute(group_id, "data", FLOAT, yields, nval);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
193 194

  /* Read the labels */
Loic Hausammann's avatar
Loic Hausammann committed
195 196
  char *labels = malloc(nval * GEAR_LABELS_SIZE);
  io_read_string_array_attribute(group_id, "elts", labels, nval,
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
197 198 199 200 201 202 203
                                 GEAR_LABELS_SIZE);

  /* Save the yields */
  /* Loop over the elements in sm */
  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
    int found = 0;
    /* Loop over SNIa yields labels */
Loic Hausammann's avatar
Loic Hausammann committed
204
    for (size_t j = 0; j < nval; j++) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
205 206 207 208 209 210 211 212 213 214 215
      const char *s1 = labels + j * GEAR_LABELS_SIZE;
      const char *s2 = stellar_evolution_get_element_name(sm, i);
      if (strcmp(s1, s2) == 0) {
        found = 1;
        snia->yields[i] = yields[j];
        break;
      }
    }

    /* Check if found an element */
    if (!found) {
Loic Hausammann's avatar
Loic Hausammann committed
216
      error("Cannot find element '%s' in SNIa yields",
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
217 218 219 220 221 222
            stellar_evolution_get_element_name(sm, i));
    }
  }

  /* Cleanup everything */
  free(yields);
Loic Hausammann's avatar
Loic Hausammann committed
223
  free(labels);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
224 225 226 227 228 229 230 231
  h5_close_group(file_id, group_id);
};

/**
 * @brief Initialize the companion structure in the #supernovae_ia.
 */
void supernovae_ia_init_companion(struct supernovae_ia *snia) {

Loic Hausammann's avatar
Loic Hausammann committed
232 233
  /* Get the exponent for the integral */
  const float exp = snia->companion_exponent;
Loic Hausammann's avatar
Loic Hausammann committed
234

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
235
  for (int i = 0; i < GEAR_NUMBER_TYPE_OF_COMPANION; i++) {
Loic Hausammann's avatar
Loic Hausammann committed
236

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
237
    /* Compute the integral */
Loic Hausammann's avatar
Loic Hausammann committed
238 239 240
    float integral = pow(snia->companion[i].mass_max, exp + 1.);
    integral -= pow(snia->companion[i].mass_min, exp + 1.);
    integral /= exp + 1.;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
241 242

    /* Update the coefficient for a normalization to 1 of the IMF */
Loic Hausammann's avatar
Loic Hausammann committed
243
    snia->companion[i].coef /= exp * integral;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
244 245 246 247 248 249 250 251
  }
}

/**
 * @brief Reads the supernovae Ia parameters from the tables.
 *
 * @param snia The #supernovae_ia model.
 * @param params The simulation parameters.
Loic Hausammann's avatar
Loic Hausammann committed
252
 * @param filename The filename of the chemistry table.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
253 254
 */
void supernovae_ia_read_from_tables(struct supernovae_ia *snia,
Loic Hausammann's avatar
Loic Hausammann committed
255 256
                                    struct swift_params *params,
                                    const char *filename) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
257 258 259 260

  hid_t file_id, group_id;

  /* Open IMF group */
261
  h5_open_group(filename, "Data/SNIa", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295

  /* Read the exponent of the IMF for companion */
  io_read_attribute(group_id, "a", FLOAT, &snia->companion_exponent);

  /* Read the minimal mass for a white dwarf progenitor */
  io_read_attribute(group_id, "Mpl", FLOAT, &snia->mass_min_progenitor);

  /* Read the maximal mass for a white dwarf */
  io_read_attribute(group_id, "Mpu", FLOAT, &snia->mass_max_progenitor);

  /* Read the maximal mass of a red giant companion */
  io_read_attribute(group_id, "Mdu1", FLOAT, &snia->companion[0].mass_max);

  /* Read the minimal mass of a red giant companion */
  io_read_attribute(group_id, "Mdl1", FLOAT, &snia->companion[0].mass_min);

  /* Read the coefficient of the main sequence companion */
  io_read_attribute(group_id, "bb1", FLOAT, &snia->companion[0].coef);

  /* Read the maximal mass of a main sequence companion */
  io_read_attribute(group_id, "Mdu2", FLOAT, &snia->companion[1].mass_max);

  /* Read the minimal mass of a main sequence companion */
  io_read_attribute(group_id, "Mdl2", FLOAT, &snia->companion[1].mass_min);

  /* Read the coefficient of the main sequence companion */
  io_read_attribute(group_id, "bb2", FLOAT, &snia->companion[1].coef);

  /* Cleanup everything */
  h5_close_group(file_id, group_id);

  /* Read the white dwarf mass */

  /* Open IMF group */
296
  h5_open_group(filename, "Data", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320

  /* Read the white dwarf mass */
  io_read_attribute(group_id, "MeanWDMass", FLOAT, &snia->mass_white_dwarf);

  /* Cleanup everything */
  h5_close_group(file_id, group_id);
}

/**
 * @brief Initialize the #supernovae_ia structure.
 *
 * @param snia The #supernovae_ia model.
 * @param phys_const The #phys_const.
 * @param us The #unit_system.
 * @param params The simulation parameters.
 * @param sm The #stellar_model.
 */
void supernovae_ia_init(struct supernovae_ia *snia,
                        const struct phys_const *phys_const,
                        const struct unit_system *us,
                        struct swift_params *params,
                        const struct stellar_model *sm) {

  /* Read the parameters from the tables */
Loic Hausammann's avatar
Loic Hausammann committed
321
  supernovae_ia_read_from_tables(snia, params, sm->yields_table);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
322 323

  /* Read the yields */
Loic Hausammann's avatar
Loic Hausammann committed
324
  supernovae_ia_read_yields(snia, params, sm, sm->yields_table);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375

  /* Get the IMF parameters */
  snia->progenitor_exponent = initial_mass_function_get_exponent(
      &sm->imf, snia->mass_min_progenitor, snia->mass_max_progenitor);
  snia->progenitor_coef_exp = initial_mass_function_get_coefficient(
      &sm->imf, snia->mass_min_progenitor, snia->mass_max_progenitor);
  snia->progenitor_coef_exp /= snia->progenitor_exponent;

  /* Compute the normalization coefficients of the companion IMF */
  supernovae_ia_init_companion(snia);
}

/**
 * @brief Write a supernovae_ia struct to the given FILE as a stream of bytes.
 *
 * Here we are only writing the arrays, everything else has been copied in the
 * feedback.
 *
 * @param snia the struct
 * @param stream the file stream
 * @param sm The #stellar_model.
 */
void supernovae_ia_dump(const struct supernovae_ia *snia, FILE *stream,
                        const struct stellar_model *sm) {

  /* Nothing to do here */
}

/**
 * @brief Restore a supernovae_ia struct from the given FILE as a stream of
 * bytes.
 *
 * Here we are only writing the arrays, everything else has been copied in the
 * feedback.
 *
 * @param snia the struct
 * @param stream the file stream
 * @param sm The #stellar_model.
 */
void supernovae_ia_restore(struct supernovae_ia *snia, FILE *stream,
                           const struct stellar_model *sm) {

  /* Nothing to do here */
}

/**
 * @brief Clean the allocated memory.
 *
 * @param snia the #supernovae_ia.
 */
void supernovae_ia_clean(struct supernovae_ia *snia) {}