Most of the cooling rates seem to be produced correctly. Need work on e.g. iron at low temperatures

examples/CoolingBox/coolingBox.yml

@@ -98,3 +98,6 @@ EAGLEChemistry:
   CalciumOverSilicon:      0.0941736
   SulphurOverSilicon:      0.6054160
 
+GearChemistry:
+  InitialMetallicity: 0.01295
+

src/chemistry/EAGLE/chemistry.h

@@ -40,10 +40,12 @@
 /**
  * @brief Return a string containing the name of a given #chemistry_element.
  */
-__attribute__((always_inline)) INLINE static const char*
+//__attribute__((always_inline)) INLINE static const char*
+__attribute__((always_inline)) INLINE const char*
 chemistry_get_element_name(enum chemistry_element elem) {
 
-  static const char* chemistry_element_names[chemistry_element_count] = {
+  //static const char* chemistry_element_names[chemistry_element_count] = {
+  const char* chemistry_element_names[chemistry_element_count] = {
       "Hydrogen", "Helium",    "Carbon",  "Nitrogen", "Oxygen",
       "Neon",     "Magnesium", "Silicon", "Iron"};
 
@@ -172,22 +174,22 @@ static INLINE void chemistry_print_backend(
  * @param p particle struct
  * @param elem enum value of element
  */
-__attribute__((always_inline)) INLINE static double chemistry_get_number_density(const struct part* restrict p, enum chemistry_element elem, const struct phys_const* restrict internal_const) {
+__attribute__((always_inline)) INLINE static double chemistry_get_number_density(const struct part* restrict p, const struct cosmology* restrict cosmo, enum chemistry_element elem, const struct phys_const* restrict internal_const) {
   double number_density;
   int atomic_number;
   switch(elem){
-    case chemistry_element_H : atomic_number = 1;
-    case chemistry_element_He: atomic_number = 4;
-    case chemistry_element_C : atomic_number = 12;
-    case chemistry_element_N : atomic_number = 14;
-    case chemistry_element_O : atomic_number = 16;
-    case chemistry_element_Ne: atomic_number = 20;
-    case chemistry_element_Mg: atomic_number = 24;
-    case chemistry_element_Si: atomic_number = 28;
-    case chemistry_element_Fe: atomic_number = 56;
+    case chemistry_element_H : atomic_number = 1; break;
+    case chemistry_element_He: atomic_number = 4; break;
+    case chemistry_element_C : atomic_number = 12; break;
+    case chemistry_element_N : atomic_number = 14; break;
+    case chemistry_element_O : atomic_number = 16; break;
+    case chemistry_element_Ne: atomic_number = 20; break;
+    case chemistry_element_Mg: atomic_number = 24; break;
+    case chemistry_element_Si: atomic_number = 28; break;
+    case chemistry_element_Fe: atomic_number = 56; break;
   }
   double element_mass = internal_const->const_proton_mass*atomic_number;
-  number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_comoving_density(p)/element_mass;
+  number_density = p->chemistry_data.metal_mass_fraction[elem]*hydro_get_physical_density(p,cosmo)/element_mass;
 
   return number_density;
 }

src/cooling/EAGLE/cooling_struct.h

@@ -99,6 +99,8 @@ struct cooling_function_data {
   float *SolarElectronAbundance;
   char **ElementNames;
   char **SolarAbundanceNames;
+  int *ElementNamePointers;
+  int *SolarAbundanceNamePointers;
 
   /*! Constant multiplication factor for time-step criterion */
   float cooling_tstep_mult;

src/cooling/EAGLE/eagle_cool_tables.h

@@ -604,6 +604,17 @@ inline struct eagle_cooling_table eagle_readtable(char *cooling_table_path, cons
   return table;
 }
 
+inline int element_index(char *element_name, const struct cooling_function_data* restrict cooling) {
+  int i;
+
+  for (i = 0; i < cooling->N_Elements; i++)
+    if (strcmp(cooling->ElementNames[i], element_name) == 0) return i;
+
+  /* element not found */
+  return -1;
+}
+
+
 inline int get_element_index(char *table[20], int size, char *element_name) {
   int i;
 
@@ -614,19 +625,153 @@ inline int get_element_index(char *table[20], int size, char *element_name) {
   return -1;
 }
 
-inline int set_cooling_SolarAbundances(const float *element_abundance,
-                                float *cooling_element_abundance, 
-				const struct cooling_function_data* restrict cooling) {
-  int i, index;
-
-  for(i = 0; i < cooling->N_Elements; i++){
-    index = get_element_index(cooling->SolarAbundanceNames,
-                          cooling->N_SolarAbundances, cooling->ElementNames[i]);
-    if (cooling->SolarAbundances[index] != 0.0) cooling_element_abundance[i] = element_abundance[i]/cooling->SolarAbundances[index];
-    else cooling_element_abundance[i] = 0.0;
+inline void MakeNamePointers(struct cooling_function_data* cooling) {
+  int i, j, sili_index = 0;
+  //char *ElementNames;
+  char ElementNames[cooling->N_Elements][eagle_element_name_length];
+
+  /* This is ridiculous, way too many element name arrays. Needs to be changed */
+  //ElementNames = malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
+  strcpy(ElementNames[0], "Hydrogen");
+  strcpy(ElementNames[1], "Helium");
+  strcpy(ElementNames[2], "Carbon");
+  strcpy(ElementNames[3], "Nitrogen");
+  strcpy(ElementNames[4], "Oxygen");
+  strcpy(ElementNames[5], "Neon");
+  strcpy(ElementNames[6], "Magnesium");
+  strcpy(ElementNames[7], "Silicon");
+  strcpy(ElementNames[8], "Iron");
+
+  cooling->ElementNamePointers = malloc(cooling->N_Elements * sizeof(int));
+  cooling->SolarAbundanceNamePointers = malloc(cooling->N_Elements * sizeof(int));
+
+  for (i = 0; i < cooling->N_Elements; i++) {
+    if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
   }
 
-  return 0;
+  for (i = 0; i < cooling->N_Elements; i++) {
+    cooling->SolarAbundanceNamePointers[i] = -999;
+    cooling->ElementNamePointers[i] = -999;
+
+    for (j = 0; j < cooling->N_SolarAbundances; j++) {
+      if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) == 0)
+        cooling->SolarAbundanceNamePointers[i] = j;
+    }
+
+    if (strcmp(cooling->ElementNames[i], "Sulphur") == 0 ||
+        strcmp(cooling->ElementNames[i], "Calcium") ==
+            0) /* These elements are tracked! */
+      cooling->ElementNamePointers[i] = -1 * sili_index;
+    else {
+      for (j = 0; j < cooling->N_Elements; j++) {
+        if (strcmp(cooling->ElementNames[i], ElementNames[j]) == 0)
+          cooling->ElementNamePointers[i] = j;
+      }
+    }
+  }
 }
 
+
+//inline int set_cooling_SolarAbundances(const float *element_abundance,
+//                                float *cooling_element_abundance, 
+//				const struct cooling_function_data* restrict cooling) {
+//  int i, index;
+//
+//  for(i = 0; i < cooling->N_Elements; i++){
+//    index = get_element_index(cooling->SolarAbundanceNames,
+//                          cooling->N_SolarAbundances, cooling->ElementNames[i]);
+//    if (cooling->SolarAbundances[index] != 0.0) cooling_element_abundance[i] = element_abundance[i]/cooling->SolarAbundances[index];
+//    else cooling_element_abundance[i] = 0.0;
+//    printf ("eagle_cool_tables.h element, name, abundance, solar abundance, cooling abundance %d %s %.5e %.5e %.5e\n",index,cooling->ElementNames[i], element_abundance[i],cooling->SolarAbundances[index], cooling_element_abundance[i]);
+//  }
+//
+//  return 0;
+//}
+
+//inline int set_cooling_SolarAbundances(const float *element_abundance,
+//                                float *cooling_element_abundance,
+//				const struct cooling_function_data* restrict cooling) {
+//  int i, index;
+//
+//  int static Silicon_SPH_Index = -1;
+//  int static Calcium_SPH_Index = -1;
+//  int static Sulphur_SPH_Index = -1;
+//
+//  int static Silicon_CoolHeat_Index = -1;
+//  int static Calcium_CoolHeat_Index = -1;
+//  int static Sulphur_CoolHeat_Index = -1;
+//
+//  static int first_call = 0;
+//
+//  if (first_call == 0) {
+//    /* determine (inverse of) solar abundance of these elements */
+//    for (i = 0; i < cooling_N_Elements; i++) {
+//      index =
+//          get_element_index(cooling_SolarAbundanceNames,
+//                            cooling_N_SolarAbundances, cooling_ElementNames[i]);
+//
+//      if (index < 0) endrun(-12345);
+//
+//      index = SolarAbundanceNamePointers[i];
+//
+//      cooling_ElementAbundance_SOLARM1[i] = 1. / cooling_SolarAbundances[index];
+//
+//      index = ElementNamePointers[i];
+//
+//      if (index < 0 && ThisTask == 0)
+//        printf("[bg_cooling] element not found %s\n", cooling_ElementNames[i]);
+//    }
+//
+//    /* Sulphur tracks Silicon: may choose not to follow Sulphur as SPH element
+//     */
+//    /* Same is true for Calcium */
+//    /* We will assume the code tracks Silicon, and may need to scale Calcium and
+//     * Sulphur accordingly */
+//
+//    Silicon_SPH_Index = element_index("Silicon");
+//    Calcium_SPH_Index = element_index("Calcium");
+//    Sulphur_SPH_Index = element_index("Sulphur");
+//
+//    Silicon_CoolHeat_Index =
+//        get_element_index(cooling_ElementNames, cooling_N_Elements, "Silicon");
+//    Calcium_CoolHeat_Index =
+//        get_element_index(cooling_ElementNames, cooling_N_Elements, "Calcium");
+//    Sulphur_CoolHeat_Index =
+//        get_element_index(cooling_ElementNames, cooling_N_Elements, "Sulphur");
+//
+//    if (Silicon_CoolHeat_Index == -1 || Calcium_CoolHeat_Index == -1 ||
+//        Sulphur_CoolHeat_Index == -1) {
+//      if (ThisTask == 0)
+//        printf("[bg_cooling] error: did not find Si or Ca or S??\n");
+//      endrun(-1233);
+//    }
+//
+//    first_call = 1;
+//  }
+//  for (i = 0; i < cooling_N_Elements; i++) {
+//    if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
+//      /* SPH does not track Calcium: use Si abundance */
+//      if (Silicon_SPH_Index == -1)
+//        cooling_element_abundance[i] = 0.0;
+//      else
+//        cooling_element_abundance[i] =
+//            element_abundance[Silicon_SPH_Index] *
+//            cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+//    else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
+//      /* SPH does not track Sulphur: use Si abundance */
+//      if (Silicon_SPH_Index == -1)
+//        cooling_element_abundance[i] = 0.0;
+//      else
+//        cooling_element_abundance[i] =
+//            element_abundance[Silicon_SPH_Index] *
+//            cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
+//    else
+//      cooling_element_abundance[i] = element_abundance[ElementNamePointers[i]] *
+//                                     cooling_ElementAbundance_SOLARM1[i];
+//  }
+//
+//  return 0;
+//}
+
+
 #endif

tests/plots.py

+import matplotlib.pyplot as plt
+
+elements = 9
+temperature = []
+cooling_rate = [[] for i in range(elements+1)]
+length = 0
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+	data = line.split()
+	temperature.append(float(data[0]))
+	cooling_rate[0].append(-float(data[1]))
+
+file_in.close()
+
+for elem in range(elements):
+	file_in = open('cooling_output_'+str(elem)+'.dat','r')
+	for line in file_in:
+        	data = line.split()
+        	cooling_rate[elem+1].append(-float(data[0]))
+	file_in.close()
+
+
+p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'b', label = 'Carbon')
+p2, = plt.loglog(temperature, cooling_rate[2], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p3, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p4, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'c', label = 'Neon')
+p5, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p6, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'y', label = 'Silicon')
+p7, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p8, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p9, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
+plt.ylim([1e-24,1e-21])
+plt.xlim([1e4,1e8])
+plt.xlabel('Temperature (K)')
+plt.ylabel('Cooling rate (eagle units...)')
+plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9])
+plt.show()

tests/testCooling.c

@@ -30,6 +30,7 @@ int main() {
   struct unit_system us;
   struct chemistry_data chemistry_data;
   struct part p;
+  struct xpart xp;
   struct phys_const internal_const;
   struct cooling_function_data cooling;
   struct cosmology cosmo;
@@ -50,43 +51,60 @@ int main() {
   units_init(&us, params, "InternalUnitSystem");
   phys_const_init(&us, params, &internal_const);
 
-  double number_density_cgs = 0.1;
-  double temperature_cgs = 1.0e6;
-  //double power_scale = units_cgs_conversion_factor(&us,UNIT_CONV_POWER)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS);
-  double power_per_num_density_factor = units_cgs_conversion_factor(&us,UNIT_CONV_POWER)*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3)/number_density_cgs;
-
-  double gamma = 5.0/3.0;
-
-  double number_density = number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
-  p.rho = number_density*(1.0/0.6*internal_const.const_proton_mass);
-  
-  for (int i = 0; i < 9; i++) p.chemistry_data.metal_mass_fraction[i] = 0.0;
-  p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen] = 0.752;
-  p.chemistry_data.metal_mass_fraction[EAGLE_Helium] = 0.248;
-
-  double temperature = temperature_cgs/units_cgs_conversion_factor(&us,UNIT_CONV_TEMPERATURE);
-  double pressure = number_density*internal_const.const_boltzmann_k*temperature;
-  printf("non-dim number density, code number density, number density scale %.5e, %.5e, %.5e\n",number_density, chemistry_get_number_density(&p,chemistry_element_H,&internal_const), pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3));
-  printf("number density, boltzmann constant, temperature: %.5e, %.5e, %.5e\n",number_density_cgs, internal_const.const_boltzmann_k, temperature);
-  printf("proton mass, boltzmann constant, %.5e, %.5e\n", internal_const.const_proton_mass, internal_const.const_boltzmann_k);
-  double internal_energy = pressure/(p.rho*(gamma - 1.0));
-  //p.entropy = internal_energy/((gamma - 1.0)*pow(p.rho,gamma - 1.0));
-  p.entropy = pressure/(pow(p.rho,gamma));
-  printf("double check pressure, actual pressure: %.5e, %.5e\n", hydro_get_comoving_pressure(&p), pressure);
-  printf("temperature, pressure, internal energy, entropy: %.5e, %.5e, %.5e, %.5e, %.5e\n", temperature,pressure,internal_energy,p.entropy,internal_const.const_boltzmann_k);
+  double hydrogen_number_density_cgs = 1e-4;
+  double u, hydrogen_number_density, pressure, gamma, cooling_du_dt, temperature_cgs;
+  int n_t_i = 2000;
+
+  char *output_filename = "cooling_output.dat";
+  FILE *output_file = fopen(output_filename, "w");
+  if (output_file == NULL)
+  {
+      printf("Error opening file!\n");
+      exit(1);
+  }
+  char *output_filename2 = "temperature_output.dat";
+  FILE *output_file2 = fopen(output_filename2, "w");
+  if (output_file2 == NULL)
+  {
+      printf("Error opening file!\n");
+      exit(1);
+  }
 
   cosmology_init(params, &us, &internal_const, &cosmo);
   cosmology_print(&cosmo);
-  printf("testCooling.c redshift %.5e\n", cosmo.z);
 
   cooling_init(params, &us, &internal_const, &cooling);
   cooling_print(&cooling);
+  //for(int j = 0; j < cooling.N_Redshifts; j++) printf("redshift %.5e\n",cooling.Redshifts[j]);
 
   chemistry_init(params, &us, &internal_const, &chemistry_data);
+  chemistry_first_init_part(&p,&xp,&chemistry_data);
   chemistry_print(&chemistry_data);
-
-  double cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const)*power_per_num_density_factor;
-  printf("cooling rate: %.5e\n",cooling_du_dt);
+  for (int k = 0; k < cooling.N_Elements; k++){
+    printf("element, abundances %d, %.5e\n", k, p.chemistry_data.metal_mass_fraction[k]);
+  }
+
+  for(int t_i = 0; t_i < n_t_i; t_i++){
+    //for (int element = 0; element < cooling.N_Elements; element++){
+    //  printf("element %d abundance %.5e particle abundance %.5e\n",element,chemistry_data.initial_metal_mass_fraction[element],p.chemistry_data.metal_mass_fraction[element]);
+    //}
+    
+    u = 1.0*pow(10.0,11 + t_i*6.0/n_t_i)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
+    pressure = u*p.rho*(gamma -1.0);
+    hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
+    p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+    p.entropy = pressure/(pow(p.rho,gamma));
+
+    gamma = 5.0/3.0;
+
+    cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const);
+    temperature_cgs = eagle_convert_u_to_temp(&p,&cooling,&cosmo,&internal_const);
+    //fprintf(output_file,"%.5e %.5e\n",u*units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS),cooling_du_dt);
+    fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
+    fprintf(output_file2,"%.5e %.5e\n",u*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)), temperature_cgs);
+  }
+  fclose(output_file);
+  fclose(output_file2);
 
   free(params);