supernovae_ii.c 13.8 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
/*******************************************************************************
 * 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_ii.h"

/* Local headers */
24
#include "engine.h"
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
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
#include "hdf5_functions.h"
#include "interpolation.h"
#include "stellar_evolution.h"
#include "stellar_evolution_struct.h"

/**
 * @brief Print the supernovae II model.
 *
 * @param snii The #supernovae_ii.
 */
void supernovae_ii_print(const struct supernovae_ii *snii) {

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

  message("Mass range for SNII = [%g, %g]", snii->mass_min, snii->mass_max);
}

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

  if (m_high < snii->mass_min || m_low > snii->mass_max) return 0;

  return 1;
}

/**
 * @brief Compute the number of supernovae II per unit of mass (equation 3.47 in
 * Poirier 2004).
 *
 * @param snii The #supernovae_ii model.
 * @param m1 The lower mass limit.
 * @param m2 The upper mass limit.
 *
 * @return The number of supernovae II per unit of mass.
 */
Loic Hausammann's avatar
Loic Hausammann committed
72 73
float supernovae_ii_get_number_per_unit_mass(const struct supernovae_ii *snii,
                                             float m1, float m2) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
#ifdef SWIFT_DEBUG_CHECKS
  if (m1 > m2) error("Mass 1 larger than mass 2 %g > %g.", m1, m2);
#endif

  /* Can we explode SNII? */
  if (!supernovae_ii_can_explode(snii, m1, m2)) {
    return 0.;
  }

  const float mass_min = max(m1, snii->mass_min);
  const float mass_max = min(m2, snii->mass_max);

  const float pow_mass =
      pow(mass_max, snii->exponent) - pow(mass_min, snii->exponent);

  return snii->coef_exp * pow_mass;
};

/**
 * @brief Get the SNII yields per mass (Poirier version).
 *
 * @param snii The #supernovae_ii model.
Loic Hausammann's avatar
Loic Hausammann committed
96 97
 * @param log_m1 The lower mass in log.
 * @param log_m2 The upper mass in log.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
98 99
 * @param yields The elements ejected (needs to be allocated).
 */
Loic Hausammann's avatar
Loic Hausammann committed
100 101 102
void supernovae_ii_get_yields_from_integral(const struct supernovae_ii *snii,
                                            float log_m1, float log_m2,
                                            float *yields) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
103 104

  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
105 106
    float yields_1 = interpolate_1d(&snii->integrated.yields[i], log_m1);
    float yields_2 = interpolate_1d(&snii->integrated.yields[i], log_m2);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
107 108 109 110 111

    yields[i] = yields_2 - yields_1;
  }
};

112 113 114 115 116 117 118
/**
 * @brief Get the SNII yields per mass.
 *
 * @param snii The #supernovae_ii model.
 * @param log_m The mass in log.
 * @param yields The elements ejected (needs to be allocated).
 */
Loic Hausammann's avatar
Loic Hausammann committed
119 120
void supernovae_ii_get_yields_from_raw(const struct supernovae_ii *snii,
                                       float log_m, float *yields) {
121 122 123 124 125 126

  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
    yields[i] = interpolate_1d(&snii->raw.yields[i], log_m);
  }
};

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
127
/**
Loic Hausammann's avatar
Loic Hausammann committed
128
 * @brief Get the ejected mass (non processed) per mass unit.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
129 130
 *
 * @param snii The #supernovae_ii model.
Loic Hausammann's avatar
Loic Hausammann committed
131 132
 * @param log_m1 The lower mass in log.
 * @param log_m2 The upper mass in log.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
133
 *
Loic Hausammann's avatar
Loic Hausammann committed
134
 * @return mass_ejected_processed The mass of non processsed elements.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
135
 */
Loic Hausammann's avatar
Loic Hausammann committed
136 137
float supernovae_ii_get_ejected_mass_fraction_non_processed_from_integral(
    const struct supernovae_ii *snii, float log_m1, float log_m2) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
138

Loic Hausammann's avatar
Loic Hausammann committed
139 140 141 142
  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);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
143 144 145 146

  return mass_ejected_2 - mass_ejected_1;
};

147
/**
Loic Hausammann's avatar
Loic Hausammann committed
148
 * @brief Get the ejected mass (non processed) per mass unit.
149 150 151 152
 *
 * @param snii The #supernovae_ii model.
 * @param log_m The mass in log.
 *
Loic Hausammann's avatar
Loic Hausammann committed
153
 * @return The mass of non processsed elements.
154
 */
Loic Hausammann's avatar
Loic Hausammann committed
155 156
float supernovae_ii_get_ejected_mass_fraction_non_processed_from_raw(
    const struct supernovae_ii *snii, float log_m) {
157

Loic Hausammann's avatar
Loic Hausammann committed
158
  return interpolate_1d(&snii->raw.ejected_mass_non_processed, log_m);
159 160
};

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
161 162 163 164 165 166 167 168 169
/**
 * @brief Get the ejected mass (processed) per mass.
 *
 * @param snii The #supernovae_ii model.
 * @param log_m1 The lower mass in log.
 * @param log_m2 The upper mass in log.
 *
 * @return mass_ejected The mass of non processsed elements.
 */
170
float supernovae_ii_get_ejected_mass_fraction_processed_from_integral(
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
171 172 173
    const struct supernovae_ii *snii, float log_m1, float log_m2) {

  float mass_ejected_1 =
174
      interpolate_1d(&snii->integrated.ejected_mass_processed, log_m1);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
175
  float mass_ejected_2 =
176
      interpolate_1d(&snii->integrated.ejected_mass_processed, log_m2);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
177 178 179 180

  return mass_ejected_2 - mass_ejected_1;
};

181 182 183 184 185 186 187 188 189 190 191 192 193 194
/**
 * @brief Get the ejected mass (processed) per mass.
 *
 * @param snii The #supernovae_ii model.
 * @param log_m The mass in log.
 *
 * @return mass_ejected The mass of non processsed elements.
 */
float supernovae_ii_get_ejected_mass_fraction_processed_from_raw(
    const struct supernovae_ii *snii, float log_m) {

  return interpolate_1d(&snii->raw.ejected_mass_processed, log_m);
};

Loic Hausammann's avatar
Gear  
Loic Hausammann committed
195 196 197 198
/**
 * @brief Read an array of SNII yields from the table.
 *
 * @param snii The #supernovae_ii model.
199 200 201 202 203 204
 * @param interp_raw Interpolation data to initialize (raw).
 * @param interp_int Interpolation data to initialize (integrated).
 * @param sm * The #stellar_model.
 * @param group_id The HDF5 group id where to read from.
 * @param hdf5_dataset_name The dataset name to read.
 * @param previous_count Number of element in the previous array read.
Loic Hausammann's avatar
Loic Hausammann committed
205 206
 * @param interpolation_size Number of element to keep in the interpolation
 * data.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
207 208
 */
void supernovae_ii_read_yields_array(
209
    struct supernovae_ii *snii, struct interpolation_1d *interp_raw,
210 211
    struct interpolation_1d *interp_int, const struct stellar_model *sm,
    hid_t group_id, const char *hdf5_dataset_name, hsize_t *previous_count,
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247
    int interpolation_size) {

  /* Now let's get the number of elements */
  /* Open attribute */
  const hid_t h_dataset = H5Dopen(group_id, hdf5_dataset_name, H5P_DEFAULT);
  if (h_dataset < 0)
    error("Error while opening attribute '%s'", hdf5_dataset_name);

  /* Get the number of elements */
  hsize_t count = io_get_number_element_in_dataset(h_dataset);

  /* Check that all the arrays have the same size */
  if (*previous_count != 0 && count != *previous_count) {
    error("The code is not able to deal with yields arrays of different size");
  }
  *previous_count = count;

  /* Read the minimal mass (in log) */
  float log_mass_min = 0;
  io_read_attribute(h_dataset, "min", FLOAT, &log_mass_min);

  /* Read the step size (log step) */
  float step_size = 0;
  io_read_attribute(h_dataset, "step", FLOAT, &step_size);

  /* Close the attribute */
  H5Dclose(h_dataset);

  /* Allocate the memory */
  float *data = (float *)malloc(sizeof(float) * count);
  if (data == NULL)
    error("Failed to allocate the SNII yields for %s.", hdf5_dataset_name);

  /* Read the dataset */
  io_read_array_dataset(group_id, hdf5_dataset_name, FLOAT, data, count);

248 249 250
  /* Initialize the raw interpolation */
  interpolate_1d_init(interp_raw, log10(snii->mass_min), log10(snii->mass_max),
                      interpolation_size, log_mass_min, step_size, count, data,
Loic Hausammann's avatar
Loic Hausammann committed
251
                      boundary_condition_zero);
252

Loic Hausammann's avatar
Loic Hausammann committed
253 254
  initial_mass_function_integrate(&sm->imf, data, count, log_mass_min,
                                  step_size);
Loic Hausammann's avatar
Loic Hausammann committed
255 256
  // TODO: decrease count in order to keep the same distance between points

257 258
  /* Initialize the integrated interpolation */
  interpolate_1d_init(interp_int, log10(snii->mass_min), log10(snii->mass_max),
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
259
                      interpolation_size, log_mass_min, step_size, count, data,
Loic Hausammann's avatar
Loic Hausammann committed
260
                      boundary_condition_const);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
261 262 263 264 265 266 267 268 269 270 271 272 273

  /* Cleanup the memory */
  free(data);
}

/**
 * @brief Read the SNII yields from the table.
 *
 * The tables are in internal units at the end of this function.
 *
 * @param snii The #supernovae_ii model.
 * @param params The simulation parameters.
 * @param sm The #stellar_model.
274
 * @param restart Are we restarting the simulation? (Is params NULL?)
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
275 276 277
 */
void supernovae_ii_read_yields(struct supernovae_ii *snii,
                               struct swift_params *params,
278 279
                               const struct stellar_model *sm,
                               const int restart) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
280 281 282 283 284

  hid_t file_id, group_id;

  hsize_t previous_count = 0;

285 286 287 288
  if (!restart) {
    snii->interpolation_size = parser_get_opt_param_int(
        params, "GEARSupernovaeII:interpolation_size", 200);
  }
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
289 290

  /* Open IMF group */
291
  h5_open_group(sm->yields_table, "Data/SNII", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
292 293 294 295 296 297 298 299

  /* Do all the elements */
  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {

    /* Get the element name */
    const char *name = stellar_evolution_get_element_name(sm, i);

    /* Read the array */
Loic Hausammann's avatar
Loic Hausammann committed
300
    supernovae_ii_read_yields_array(
301 302
        snii, &snii->raw.yields[i], &snii->integrated.yields[i], sm, group_id,
        name, &previous_count, snii->interpolation_size);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
303 304 305
  }

  /* Read the mass ejected */
Loic Hausammann's avatar
Loic Hausammann committed
306
  supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass_processed,
307 308 309
                                  &snii->integrated.ejected_mass_processed, sm,
                                  group_id, "Ej", &previous_count,
                                  snii->interpolation_size);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
310 311

  /* Read the mass ejected of non processed gas */
Loic Hausammann's avatar
Loic Hausammann committed
312
  supernovae_ii_read_yields_array(snii, &snii->raw.ejected_mass_non_processed,
Loic Hausammann's avatar
Loic Hausammann committed
313
                                  &snii->integrated.ejected_mass_non_processed,
314 315
                                  sm, group_id, "Ejnp", &previous_count,
                                  snii->interpolation_size);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
316 317 318 319 320 321 322 323 324 325

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

/**
 * @brief Reads the supernovae II parameters from tables.
 *
 * @param snii The #supernovae_ii model.
 * @param params The simulation parameters.
Loic Hausammann's avatar
Loic Hausammann committed
326
 * @param filename The filename of the chemistry table.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
327 328
 */
void supernovae_ii_read_from_tables(struct supernovae_ii *snii,
Loic Hausammann's avatar
Loic Hausammann committed
329 330
                                    struct swift_params *params,
                                    const char *filename) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
331 332 333 334

  hid_t file_id, group_id;

  /* Open IMF group */
335
  h5_open_group(filename, "Data/SNII", &file_id, &group_id);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353

  /* Read the minimal mass of a supernovae */
  io_read_attribute(group_id, "Mmin", FLOAT, &snii->mass_min);

  /* Read the maximal mass of a supernovae */
  io_read_attribute(group_id, "Mmax", FLOAT, &snii->mass_max);

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

/**
 * @brief Initialize the #supernovae_ii structure.
 *
 * @param snii The #supernovae_ii model.
 * @param params The simulation parameters.
 * @param sm The #stellar_model.
 */
354
void supernovae_ii_init(struct supernovae_ii *snii, struct swift_params *params,
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
355 356 357
                        const struct stellar_model *sm) {

  /* Read the parameters from the tables */
Loic Hausammann's avatar
Loic Hausammann committed
358
  supernovae_ii_read_from_tables(snii, params, sm->yields_table);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
359

360 361
  /* Read the supernovae yields */
  supernovae_ii_read_yields(snii, params, sm, /* restart */ 0);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378

  /* Get the IMF parameters */
  snii->exponent = initial_mass_function_get_exponent(&sm->imf, snii->mass_min,
                                                      snii->mass_max);
  snii->coef_exp = initial_mass_function_get_coefficient(
      &sm->imf, snii->mass_min, snii->mass_max);
  snii->coef_exp /= snii->exponent;
}

/**
 * @brief Write a supernovae_ii 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 snii the struct
 * @param stream the file stream
Loic Hausammann's avatar
Loic Hausammann committed
379
 * @param sm The #stellar_model.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
380 381
 */
void supernovae_ii_dump(const struct supernovae_ii *snii, FILE *stream,
382
                        const struct stellar_model *sm) {}
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
383 384 385 386 387 388 389 390 391 392

/**
 * @brief Restore a supernovae_ii 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 snii the struct
 * @param stream the file stream
Loic Hausammann's avatar
Loic Hausammann committed
393
 * @param sm The #stellar_model.
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
394 395
 */
void supernovae_ii_restore(struct supernovae_ii *snii, FILE *stream,
396
                           const struct stellar_model *sm) {
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
397

398
  /* Read the supernovae yields (and apply the units) */
399
  supernovae_ii_read_yields(snii, NULL, sm, /* restart */ 1);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
400

401 402 403 404
  /* Get the IMF parameters */
  snii->exponent = initial_mass_function_get_exponent(&sm->imf, snii->mass_min,
                                                      snii->mass_max);
  snii->coef_exp = initial_mass_function_get_coefficient(
405
      &sm->imf, snii->mass_min, snii->mass_max);
406
  snii->coef_exp /= snii->exponent;
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
407 408 409 410 411 412 413 414 415 416
}

/**
 * @brief Clean the allocated memory.
 *
 * @param snii the #supernovae_ii.
 */
void supernovae_ii_clean(struct supernovae_ii *snii) {

  for (int i = 0; i < GEAR_CHEMISTRY_ELEMENT_COUNT; i++) {
417 418
    interpolate_1d_free(&snii->integrated.yields[i]);
    interpolate_1d_free(&snii->raw.yields[i]);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
419 420
  }

421 422
  interpolate_1d_free(&snii->integrated.ejected_mass_processed);
  interpolate_1d_free(&snii->raw.ejected_mass_processed);
Loic Hausammann's avatar
Loic Hausammann committed
423 424
  interpolate_1d_free(&snii->integrated.ejected_mass_non_processed);
  interpolate_1d_free(&snii->raw.ejected_mass_non_processed);
Loic Hausammann's avatar
Gear  
Loic Hausammann committed
425
}