testFeedback.c 11.6 KB
Newer Older
Alexei Borissov's avatar
Alexei Borissov committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
 *
 * 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 "swift.h"

22
#if defined(STARS_EAGLE) && defined(FEEDBACK_EAGLE)
23
/**
24 25 26 27
 * @brief compute the relative error between two floats
 *
 * @param a, b numbers between which to compute the relative difference
 */
28
float relative_error(float a, float b) { return fabs((a - b) / b); }
29 30 31

/**
 * @brief Test function to check feedback is working correctly. Produces a table
32 33 34
 * of values for the cumulative amount of mass enrichment up to a given time
 * from the birth of the star (normalized by the initial stellar mass). Compares
 * these results to those produced by an analogous test in EAGLE.
35
 */
Alexei Borissov's avatar
Alexei Borissov committed
36
int main(int argc, char *argv[]) {
37

Alexei Borissov's avatar
Alexei Borissov committed
38 39 40 41 42 43 44 45 46 47 48
  /* Declare relevant structs */
  struct swift_params *params = malloc(sizeof(struct swift_params));
  struct unit_system us;
  struct chemistry_global_data chem_data;
  struct part p;
  struct xpart xp;
  struct spart sp;
  struct phys_const phys_const;
  struct cosmology cosmo;
  struct hydro_props hydro_properties;
  struct stars_props stars_properties;
49
  struct feedback_props feedback_properties;
Alexei Borissov's avatar
Alexei Borissov committed
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73
  char *parametersFileName = "./testFeedback.yml";

  /* Read the parameter file */
  if (params == NULL) error("Error allocating memory for the parameter file.");
  message("Reading runtime parameters from file '%s'", parametersFileName);
  parser_read_file(parametersFileName, params);

  /* Init units */
  units_init_from_params(&us, params, "InternalUnitSystem");
  phys_const_init(&us, params, &phys_const);

  /* Init chemistry */
  chemistry_init(params, &us, &phys_const, &chem_data);
  chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
  chemistry_print(&chem_data);

  /* Init cosmology */
  cosmology_init(params, &us, &phys_const, &cosmo);
  cosmology_print(&cosmo);

  /* Init hydro properties */
  hydro_props_init(&hydro_properties, &phys_const, &us, params);

  /* Init star properties */
74 75
  stars_props_init(&stars_properties, &phys_const, &us, params,
                   &hydro_properties, &cosmo);
Alexei Borissov's avatar
Alexei Borissov committed
76

77 78 79
  /* Init star properties */
  feedback_props_init(&feedback_properties, &phys_const, &us, params,
                      &hydro_properties, &cosmo);
Alexei Borissov's avatar
Alexei Borissov committed
80 81

  /* Init spart */
82 83
  stars_first_init_spart(&sp, &stars_properties, /*with_cosmology=*/0, cosmo.a,
                         cosmo.time);
Alexei Borissov's avatar
Alexei Borissov committed
84

85 86 87
  /* Define an initial stellar mass. (for use when calling the feedback
   * functions, the results are presented per initial stellar mass, so the
   * actual value does not matter. */
Alexei Borissov's avatar
Alexei Borissov committed
88 89
  sp.mass_init = 4.706273e-5;

90 91 92 93 94 95 96
  /* Set metal mass fractions */
  for (int i = 0; i < chemistry_element_count; i++)
    sp.chemistry_data.metal_mass_fraction[i] = 0.f;
  sp.chemistry_data.metal_mass_fraction[0] = 0.752;
  sp.chemistry_data.metal_mass_fraction[1] = 0.248;
  sp.chemistry_data.metal_mass_fraction_total = 0.01;

97
  /* Define how long to run for and what interval should be between consecutive
Alexei Borissov's avatar
Alexei Borissov committed
98
   * times for which feedback is calculated for */
Alexei Borissov's avatar
Alexei Borissov committed
99
  float Gyr_to_s = 3.154e16;
100 101 102
  float dt = 0.1 * Gyr_to_s / units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
  float max_age =
      13.f * Gyr_to_s / units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
103

Alexei Borissov's avatar
Alexei Borissov committed
104
  /* Zero feedback quantities */
105
  for (int i = 0; i < chemistry_element_count; i++)
106 107 108 109 110 111 112 113 114 115
    sp.feedback_data.to_distribute.metal_mass[i] = 0.f;
  sp.feedback_data.to_distribute.metal_mass_from_SNIa = 0.f;
  sp.feedback_data.to_distribute.metal_mass_from_SNII = 0.f;
  sp.feedback_data.to_distribute.metal_mass_from_AGB = 0.f;
  sp.feedback_data.to_distribute.mass_from_AGB = 0.f;
  sp.feedback_data.to_distribute.mass_from_SNII = 0.f;
  sp.feedback_data.to_distribute.mass_from_SNIa = 0.f;
  sp.feedback_data.to_distribute.Fe_mass_from_SNIa = 0.f;
  sp.feedback_data.to_distribute.total_metal_mass = 0.f;
  sp.feedback_data.to_distribute.mass = 0.f;
116

117 118
  /* Open EAGLE test file for reading  */
  FILE *EAGLE_test;
119 120
  const char EAGLE_fname[75] =
      "/cosma/home/dp004/dc-bori1/Eagle/data1/z_0.01/StellarEvolutionTotal.txt";
121 122 123
  if (!(EAGLE_test = fopen(EAGLE_fname, "r"))) {
    error("error in opening file '%s'\n", EAGLE_fname);
  }
Alexei Borissov's avatar
Alexei Borissov committed
124

125 126 127 128 129 130 131 132 133 134
  /* Declare constants necessary for reading EAGLE data */
  char *line = NULL;
  size_t len = 0;
  const int n_fields = 19;
  float tol = 1e-5;

  /* Declare array to store one line of data from EAGLE test */
  float eagle_data[n_fields];

  /* read first line */
135 136
  if (getline(&line, &len, EAGLE_test) == -1)
    error("failed to read first line of EAGLE test file");
137 138

  /* Open file for writing SWIFT feedback test output */
Alexei Borissov's avatar
Alexei Borissov committed
139 140 141 142
  FILE *Total_output;
  const char Total_fname[25] = "test_feedback_total.txt";
  if (!(Total_output = fopen(Total_fname, "w"))) {
    error("error in opening file '%s'\n", Total_fname);
Alexei Borissov's avatar
Alexei Borissov committed
143
  }
Alexei Borissov's avatar
Alexei Borissov committed
144
  fprintf(Total_output,
Alexei Borissov's avatar
Alexei Borissov committed
145
          "# time[Gyr] | total mass | metal mass: total | H | He | C | N  | O  "
146 147
          "| Ne | Mg | Si | Fe | per solar mass (m,z)_AGB (m,z)_SNII "
          "(m,z,M_fe)_SNIa \n");
Alexei Borissov's avatar
Alexei Borissov committed
148

Alexei Borissov's avatar
Alexei Borissov committed
149
  /* Loop over times for which to calculate feedback */
Alexei Borissov's avatar
Alexei Borissov committed
150
  for (float age = 0; age <= max_age; age += dt) {
Alexei Borissov's avatar
Alexei Borissov committed
151 152

    /* Compute feedback */
153
    compute_stellar_evolution(&feedback_properties, &cosmo, &sp, &us, age, dt);
Alexei Borissov's avatar
Alexei Borissov committed
154

Alexei Borissov's avatar
Alexei Borissov committed
155
    /* Print computed values to file */
Alexei Borissov's avatar
Alexei Borissov committed
156 157
    float age_Gyr =
        age * units_cgs_conversion_factor(&us, UNIT_CONV_TIME) / Gyr_to_s;
Alexei Borissov's avatar
Alexei Borissov committed
158
    fprintf(Total_output, "%f %e %e ", age_Gyr,
159 160
            sp.feedback_data.to_distribute.mass / sp.mass_init,
            sp.feedback_data.to_distribute.total_metal_mass / sp.mass_init);
Alexei Borissov's avatar
Alexei Borissov committed
161
    for (int i = 0; i < chemistry_element_count; i++)
162
      fprintf(Total_output, "%e ",
163
              sp.feedback_data.to_distribute.metal_mass[i] / sp.mass_init);
164
    fprintf(Total_output, " %e %e %e %e %e %e %e",
165 166 167 168 169 170 171
            sp.feedback_data.to_distribute.mass_from_AGB / sp.mass_init,
            sp.feedback_data.to_distribute.metal_mass_from_AGB / sp.mass_init,
            sp.feedback_data.to_distribute.mass_from_SNII / sp.mass_init,
            sp.feedback_data.to_distribute.metal_mass_from_SNII / sp.mass_init,
            sp.feedback_data.to_distribute.mass_from_SNIa / sp.mass_init,
            sp.feedback_data.to_distribute.metal_mass_from_SNIa / sp.mass_init,
            sp.feedback_data.to_distribute.Fe_mass_from_SNIa / sp.mass_init);
Alexei Borissov's avatar
Alexei Borissov committed
172
    fprintf(Total_output, "\n");
173

Alexei Borissov's avatar
Alexei Borissov committed
174
    /* Read data from EAGLE test and compare it to what we calculated */
175
    if (!feof(EAGLE_test)) {
176 177 178 179 180 181 182 183 184
      int ret = fscanf(
          EAGLE_test,
          "%e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ",
          &eagle_data[0], &eagle_data[1], &eagle_data[2], &eagle_data[3],
          &eagle_data[4], &eagle_data[5], &eagle_data[6], &eagle_data[7],
          &eagle_data[8], &eagle_data[9], &eagle_data[10], &eagle_data[11],
          &eagle_data[12], &eagle_data[13], &eagle_data[14], &eagle_data[15],
          &eagle_data[16], &eagle_data[17], &eagle_data[18]);
      if (ret == 0) error("Error reading input.");
185 186 187 188
      if (relative_error(age_Gyr, eagle_data[0]) > tol)
        error(
            "relative error in age greater than tolerance. Swift: %e Eagle %e",
            age_Gyr, eagle_data[0]);
189 190
      if (relative_error(sp.feedback_data.to_distribute.mass / sp.mass_init,
                         eagle_data[1]) > tol)
191 192 193
        error(
            "relative error in total mass greater than tolerance. Swift: %e "
            "Eagle %e",
194 195 196 197
            sp.feedback_data.to_distribute.mass / sp.mass_init, eagle_data[1]);
      if (relative_error(
              sp.feedback_data.to_distribute.total_metal_mass / sp.mass_init,
              eagle_data[2]) > tol)
198 199 200
        error(
            "relative error in total metal mass greater than tolerance. Swift: "
            "%e Eagle %e",
201 202
            sp.feedback_data.to_distribute.total_metal_mass / sp.mass_init,
            eagle_data[2]);
203
      for (int i = 0; i < chemistry_element_count; i++)
204 205 206
        if (relative_error(
                sp.feedback_data.to_distribute.metal_mass[i] / sp.mass_init,
                eagle_data[3 + i]) > tol)
207 208 209 210
          error(
              "relative error in mass released for element %d greater than "
              "tolerance. Swift: %e Eagle %e",
              i, age_Gyr, eagle_data[3 + i]);
211 212 213
      if (relative_error(
              sp.feedback_data.to_distribute.mass_from_AGB / sp.mass_init,
              eagle_data[12]) > tol)
214 215 216
        error(
            "relative error in mass from AGB greater than tolerance. Swift: %e "
            "Eagle %e",
217 218 219 220 221
            sp.feedback_data.to_distribute.mass_from_AGB / sp.mass_init,
            eagle_data[12]);
      if (relative_error(
              sp.feedback_data.to_distribute.metal_mass_from_AGB / sp.mass_init,
              eagle_data[13]) > tol)
222 223 224
        error(
            "relative error in metal mass from AGB greater than tolerance. "
            "Swift: %e Eagle %e",
225
            sp.feedback_data.to_distribute.metal_mass_from_AGB / sp.mass_init,
226
            eagle_data[13]);
227 228 229
      if (relative_error(
              sp.feedback_data.to_distribute.mass_from_SNII / sp.mass_init,
              eagle_data[14]) > tol)
230 231 232
        error(
            "relative error in mass from SNII greater than tolerance. Swift: "
            "%e Eagle %e",
233 234 235 236
            sp.feedback_data.to_distribute.mass_from_SNII / sp.mass_init,
            eagle_data[14]);
      if (relative_error(sp.feedback_data.to_distribute.metal_mass_from_SNII /
                             sp.mass_init,
237 238 239 240
                         eagle_data[15]) > tol)
        error(
            "relative error in metal mass from SNII greater than tolerance. "
            "Swift: %e Eagle %e",
241
            sp.feedback_data.to_distribute.metal_mass_from_SNII / sp.mass_init,
242
            eagle_data[15]);
243 244 245
      if (relative_error(
              sp.feedback_data.to_distribute.mass_from_SNIa / sp.mass_init,
              eagle_data[16]) > tol)
246 247 248
        error(
            "relative error in mass from SNIa greater than tolerance. Swift: "
            "%e Eagle %e",
249 250 251 252
            sp.feedback_data.to_distribute.mass_from_SNIa / sp.mass_init,
            eagle_data[16]);
      if (relative_error(sp.feedback_data.to_distribute.metal_mass_from_SNIa /
                             sp.mass_init,
253 254 255 256
                         eagle_data[17]) > tol)
        error(
            "relative error in metal mass from SNIa greater than tolerance. "
            "Swift: %e Eagle %e",
257
            sp.feedback_data.to_distribute.metal_mass_from_SNIa / sp.mass_init,
258
            eagle_data[17]);
259 260 261
      if (relative_error(
              sp.feedback_data.to_distribute.Fe_mass_from_SNIa / sp.mass_init,
              eagle_data[18]) > tol)
262 263 264
        error(
            "relative error in iron mass from SNIa greater than tolerance. "
            "Swift: %e Eagle %e",
265 266
            sp.feedback_data.to_distribute.Fe_mass_from_SNIa / sp.mass_init,
            eagle_data[18]);
267 268 269
    } else {
      error("Failed to read line of EAGLE test file data");
    }
Alexei Borissov's avatar
Alexei Borissov committed
270
  }
271

Alexei Borissov's avatar
Alexei Borissov committed
272 273
  return 0;
}
274 275 276 277

#else

/* Don't do anything if not using EAGLE stars */
278
int main(int argc, char *argv[]) { return 0; }
279 280

#endif