collision and photodissociation cooling reading in well, redshift cooling needs work

examples/CoolingBox/analytical_test.py

@@ -86,8 +86,8 @@ temp_0 = 1.0e7
 rho = rho*unit_mass/(unit_length**3)
 
 # Read snapshots
-nsnap = 40
-npart = 32768
+nsnap = 150
+npart = 4096
 u_snapshots_cgs = zeros(nsnap)
 u_part_snapshots_cgs = zeros((nsnap,npart))
 t_snapshots_cgs = zeros(nsnap)
@@ -107,7 +107,7 @@ for line in file_in:
         temperature.append(float(data[0]))
         cooling_rate.append(-float(data[1]))
 
-tfinal = 3.3*t_snapshots_cgs[nsnap-1]
+tfinal = t_snapshots_cgs[nsnap-1]
 nt = 1e4
 dt = tfinal/nt
 
@@ -121,10 +121,11 @@ for j in range(int(nt-1)):
 	t_sol[j+1] = t_sol[j] + dt
 	Lambda_net = interpol_lambda(temperature,cooling_rate,temp_sol[j])
 	#u_next = (u*m_p - Lambda_net*rho/(0.59*m_p)*dt)/m_p
-	nH = 0.75*rho/(m_p)
+	nH = 0.7*rho/(m_p)
+	#nH = 9.125e-5
 	u_next = u - Lambda_net*nH**2/rho*dt
 	#print(u_next, u, Lambda_net,rho/(0.59*m_p),dt)
-	print(Lambda_net,rho,dt,nH)
+	#print(nH**2/rho*Lambda_net, dt)
 	temp_sol[j+1] = convert_u_to_temp_sol(u_next,rho)
 	lambda_sol[j] = Lambda_net
 	u = u_next
@@ -141,9 +142,9 @@ ylabel("${\\rm{Temperature~[K]}}$")
 
 savefig("energy.png", dpi=200)
 
-p = figure()
-p1, = loglog(temp_sol,lambda_sol,linewidth = 0.5,color = 'k',label = 'analytical')
-xlabel("${\\rm{Temperature~[K]}}$")
-ylabel("${\\rm{Cooling~rate}}$")
+#p = figure()
+#p1, = loglog(temp_sol,lambda_sol,linewidth = 0.5,color = 'k',label = 'analytical')
+#xlabel("${\\rm{Temperature~[K]}}$")
+#ylabel("${\\rm{Cooling~rate}}$")
 
-savefig("cooling.png", dpi=200)
+#savefig("cooling.png", dpi=200)

examples/CoolingBox/coolingBox.yml

@@ -17,7 +17,7 @@ Cosmology:
 # Parameters governing the time integration
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
-  time_end:   0.02  # The end time of the simulation (in internal units).
+  time_end:   1.5  # The end time of the simulation (in internal units).
   dt_min:     1e-5  # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
@@ -25,7 +25,7 @@ TimeIntegration:
 Snapshots:
   basename:            coolingBox # Common part of the name of output files
   time_first:          0.         # Time of the first output (in internal units)
-  delta_time:          5e-4       # Time difference between consecutive outputs (in internal units)
+  delta_time:          1e-2       # Time difference between consecutive outputs (in internal units)
 
 # Parameters governing the conserved quantities statistics
 Statistics:
@@ -71,15 +71,15 @@ EagleCooling:
 
 EAGLEChemistry:
   InitMetallicity:         0.014
-  InitAbundance_Hydrogen:  0.706
-  InitAbundance_Helium:    0.280
-  InitAbundance_Carbon:    0.00207
-  InitAbundance_Nitrogen:  8.35e-4
-  InitAbundance_Oxygen:    0.00549
-  InitAbundance_Neon:      0.00141
-  InitAbundance_Magnesium: 5.91e-4
-  InitAbundance_Silicon:   6.83e-4
-  InitAbundance_Iron:      0.0011
+  InitAbundance_Hydrogen:  0.70649785
+  InitAbundance_Helium:    0.28055534
+  InitAbundance_Carbon:    2.0665436e-3
+  InitAbundance_Nitrogen:  8.3562563e-4
+  InitAbundance_Oxygen:    5.4926244e-3
+  InitAbundance_Neon:      1.4144605e-3
+  InitAbundance_Magnesium: 5.907064e-4
+  InitAbundance_Silicon:   6.825874e-4
+  InitAbundance_Iron:      1.1032152e-3
   CalciumOverSilicon:      0.0941736
   SulphurOverSilicon:      0.6054160
   #InitMetallicity:         0.

examples/CoolingBox/getGlass.sh

 #!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_16.hdf5

examples/CoolingBox/makeIC.py

@@ -28,7 +28,7 @@ from numpy import *
 periodic= 1           # 1 For periodic box
 boxSize = 1           # 1 kiloparsec    
 rho = 3.2e3           # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 4.5e8             # Pressure in code units (at 10^7K)
+P = 4.5e7             # Pressure in code units (at 10^6K)
 gamma = 5./3.         # Gas adiabatic index
 eta = 1.2349          # 48 ngbs with cubic spline kernel
 fileName = "coolingBox.hdf5" 
@@ -36,7 +36,7 @@ fileName = "coolingBox.hdf5"
 #---------------------------------------------------
 
 # Read id, position and h from glass
-glass = h5py.File("glassCube_32.hdf5", "r")
+glass = h5py.File("glassCube_16.hdf5", "r")
 ids = glass["/PartType0/ParticleIDs"][:]
 pos = glass["/PartType0/Coordinates"][:,:] * boxSize
 h = glass["/PartType0/SmoothingLength"][:] * boxSize

examples/CoolingBox/run.sh

@@ -2,7 +2,7 @@
 #!/bin/bash
 
 # Generate the initial conditions if they are not present.
-if [ ! -e glassCube_32.hdf5 ]
+if [ ! -e glassCube_16.hdf5 ]
 then
     echo "Fetching initial glass file for the cooling box example..."
     ./getGlass.sh
@@ -21,7 +21,7 @@ then
 fi
 
 # Run SWIFT
-../swift -s -C -t 4 coolingBox.yml
+../swift -s -C -t 8 coolingBox.yml
 
 # Check energy conservation and cooling rate
 python energy_plot.py

src/cooling/EAGLE/cooling.h

@@ -110,27 +110,6 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable
     *dx = 1;
   } else {
     *i = (int)floor(((float)x - table[0]) * dxm1);
-    //printf("Eagle cooling.h ntable, i, x, table[0], table[i], table[ntable-1], dxm1 %d %d %.5e %.5e %.5e %.5e %.5e\n", ntable, *i, x, table[0], table[*i], table[ntable-1], dxm1);
-    fflush(stdout);
-    *dx = ((float)x - table[*i]) * dxm1;
-  }
-}
-
-__attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int ntable, double x, int *i, float *dx) {
-  float dxm1;
-  const float EPS = 1.e-4;
-
-  dxm1 = (float)(ntable - 1) / (table[ntable - 1] - table[0]);
-
-  if ((float)x <= table[0] + EPS) {
-    *i = 0;
-    *dx = 0;
-  } else if ((float)x >= table[ntable - 1] - EPS) {
-    *i = ntable - 2;
-    *dx = 1;
-  } else {
-    *i = (int)floor(((float)x - table[0]) * dxm1);
-    //printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
     *dx = ((float)x - table[*i]) * dxm1;
   }
 }
@@ -185,8 +164,6 @@ __attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int
   result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
            dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
 
-  //printf("Eagle cooling.h interpol_2d dx dy table indices (0-3) nx ny %.5e %.5e %d %d %d %d %d %d\n", dx, dy, index[0], index[1], index[2], index[3], nx, ny);
-
   return result;
 }
 
@@ -243,6 +220,17 @@ __attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int
   if(index[6] >= nx*ny*nz || index[6] < 0) fprintf(stderr,"index 6 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[6],i+1,j+1,k,nx*ny*nz);
   if(index[7] >= nx*ny*nz || index[7] < 0) fprintf(stderr,"index 7 out of bounds %d, i,j,k %d, %d, %d, table size %d\n", index[7],i+1,j+1,k+1,nx*ny*nz);
 
+  //printf("index 0 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[0]], index[0],i,j,k,nx*ny*nz);
+  //printf("index 1 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[1]], index[1],i,j,k+1,nx*ny*nz);
+  //printf("index 2 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[2]], index[2],i,j+1,k,nx*ny*nz);
+  //printf("index 3 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[3]], index[3],i,j+1,k+1,nx*ny*nz);
+  //printf("index 4 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[4]], index[4],i+1,j,k,nx*ny*nz);
+  //printf("index 5 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[5]], index[5],i+1,j,k+1,nx*ny*nz);
+  //printf("index 6 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[6]], index[6],i+1,j+1,k,nx*ny*nz);
+  //printf("index 7 out of bounds %.5e %d, i,j,k %d, %d, %d, table size %d\n", table[index[7]], index[7],i+1,j+1,k+1,nx*ny*nz);
+  //fflush(stdout);
+
+
   result = (1 - dx) * (1 - dy) * (1 - dz) * table[index[0]] +
            (1 - dx) * (1 - dy) * dz * table[index[1]] +
            (1 - dx) * dy * (1 - dz) * table[index[2]] +
@@ -409,7 +397,7 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
         cooling_element_abundance[i] = 0.0;
       else{
         cooling_element_abundance[i] =
-            element_abundance[Silicon_SPH_Index] * cooling->calcium_over_silicon_ratio *
+            element_abundance[sili_index] * cooling->calcium_over_silicon_ratio *
             cooling_ElementAbundance_SOLARM1[Calcium_CoolHeat_Index];
       }
     else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index != -1)
@@ -418,7 +406,7 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
         cooling_element_abundance[i] = 0.0;
       else{
         cooling_element_abundance[i] =
-            element_abundance[Silicon_SPH_Index] * cooling->sulphur_over_silicon_ratio *
+            element_abundance[sili_index] * cooling->sulphur_over_silicon_ratio *
             cooling_ElementAbundance_SOLARM1[Sulphur_CoolHeat_Index];
       }
     else{
@@ -480,6 +468,7 @@ __attribute__((always_inline)) INLINE void get_redshift_index(float z, int *z_in
 
         break;
       }
+    //printf("Eagle cooling.h z, z_index, z_index_previous, redshifts grid elem, dz %.5e %d %d %.5e %.5e %.5e\n", z,iz,get_redshift_index_previous,*dz, cooling->Redshifts[iz], cooling->Redshifts[iz+1]);
     }
   }
 }
@@ -523,8 +512,8 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons
  * @brief calculates cooling rate
  *
  */
-__attribute__((always_inline)) INLINE static double eagle_cooling_rate(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const) {
-  double T_gam, LambdaNet = 0.0, solar_electron_abundance;
+__attribute__((always_inline)) INLINE static void eagle_metal_cooling_rate(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const, double* element_lambda) {
+  double T_gam, solar_electron_abundance;
   double n_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,internal_const)*cooling->number_density_scale; // chemistry_data
   double z = cosmo->z;
   float dz; 
@@ -546,28 +535,34 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
   get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
   get_index_1d(cooling->nH, cooling->N_nH, log10(n_h), &n_h_i, &d_n_h);
 
-  char output_filename[21];
-  FILE** output_file = malloc(9*sizeof(FILE*));
-  for (int element = 0; element < cooling->N_Elements; element++){
-    sprintf(output_filename, "%s%d%s", "cooling_output_",element,".dat");
-    output_file[element] = fopen(output_filename, "a");
-    if (output_file == NULL)
-    {   
-        printf("Error opening file!\n");
-        exit(1);
-    }
-  }
-
   /* ------------------ */
   /* Metal-free cooling */
   /* ------------------ */
 
-    LambdaNet =
-        interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, z_index, He_i, n_h_i,
-                     temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+    /* Collisional cooling */
+    //element_lambda[0] =
+    //    interpol_2d(cooling->table.collisional_cooling.H_plus_He_heating, He_i,
+    //                 temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+    //h_plus_he_electron_abundance =
+    //    interpol_2d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, He_i,
+    //                 temp_i, d_He, d_temp,cooling->N_He,cooling->N_Temp);
+    
+    /* Photodissociation */
+    element_lambda[0] =
+        interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_heating, He_i,
+                     temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
     h_plus_he_electron_abundance =
-        interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, z_index, He_i, n_h_i,
-                     temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+        interpol_3d(cooling->table.photodissociation_cooling.H_plus_He_electron_abundance, He_i,
+                     temp_i, n_h_i, d_He, d_temp, d_n_h,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+
+    /* redshift tables */
+    //printf("Eagle cooling.h z_index, dz, z %d %.5e %.5e\n", z_index, dz, cooling->Redshifts[z_index]);
+    //element_lambda[0] =
+    //    interpol_4d(cooling->table.element_cooling.H_plus_He_heating, z_index, He_i,
+    //                 temp_i, n_h_i, dz, d_He, d_temp, d_n_h, cooling->N_Redshifts,cooling->N_He,cooling->N_Temp,cooling->N_nH);
+    //h_plus_he_electron_abundance =
+    //    interpol_4d(cooling->table.element_cooling.H_plus_He_electron_abundance, z_index, He_i,
+    //                 temp_i, n_h_i, dz, d_He, d_temp, d_n_h, cooling->N_Redshifts,cooling->N_He,cooling->N_Temp,cooling->N_nH);
 
   /* ------------------ */
   /* Compton cooling    */
@@ -580,7 +575,7 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
 
     T_gam = eagle_cmb_temperature * (1 + z);
 
-    LambdaNet -= eagle_compton_rate * (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
+    element_lambda[1] = -eagle_compton_rate * (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
                  h_plus_he_electron_abundance / n_h;
   }
 
@@ -593,24 +588,78 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
 
     set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling, p);
 
+    /* Collisional cooling */
+    //solar_electron_abundance =
+    //    interpol_1d(cooling->table.collisional_cooling.electron_abundance, temp_i, d_temp); /* ne/n_h */
+
+    //for (i = 0; i < cooling->N_Elements; i++){
+    //    element_lambda[i+2] = interpol_2d(cooling->table.collisional_cooling.metal_heating, i,
+    //                    temp_i, 0.0, d_temp,cooling->N_Elements,cooling->N_Temp) *
+    //        (h_plus_he_electron_abundance / solar_electron_abundance) *
+    //        cooling_solar_abundances[i];
+    //}
+    
+    /* Photodissociation */
     solar_electron_abundance =
-        interpol_3d(cooling->table.element_cooling.electron_abundance, z_index, n_h_i, temp_i, dz,
-                    d_n_h, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); /* ne/n_h */
+        interpol_2d(cooling->table.photodissociation_cooling.electron_abundance, temp_i, n_h_i, d_temp, d_n_h, cooling->N_Temp, cooling->N_nH); /* ne/n_h */
 
-    double element_cooling_lambda;
     for (i = 0; i < cooling->N_Elements; i++){
-        element_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i, n_h_i,
-                        temp_i, dz, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
+        element_lambda[i+2] = interpol_3d(cooling->table.photodissociation_cooling.metal_heating, i,
+                        temp_i, n_h_i, 0.0, d_temp, d_n_h,cooling->N_Elements,cooling->N_Temp,cooling->N_nH) *
             (h_plus_he_electron_abundance / solar_electron_abundance) *
             cooling_solar_abundances[i];
-        LambdaNet += element_cooling_lambda;
-	fprintf(output_file[i],"%.5e\n",element_cooling_lambda);
     }
+    
+    /* redshift tables */
+    //solar_electron_abundance =
+    //    interpol_3d(cooling->table.element_cooling.electron_abundance, z_index, temp_i, n_h_i, dz, d_temp, d_n_h, cooling->N_Redshifts, cooling->N_Temp, cooling->N_nH); /* ne/n_h */
+
+    //for (i = 0; i < cooling->N_Elements; i++){
+    //    element_lambda[i+2] = interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i,
+    //                    temp_i, n_h_i, dz, 0.0, d_temp, d_n_h,cooling->N_Redshifts,cooling->N_Elements,cooling->N_Temp,cooling->N_nH) *
+    //        (h_plus_he_electron_abundance / solar_electron_abundance) *
+    //        cooling_solar_abundances[i];
+    //}
+}
 
 
-  for (i = 0; i < cooling->N_Elements; i++) fclose(output_file[i]);
+__attribute__((always_inline)) INLINE static double eagle_cooling_rate(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const) {
+  double *element_lambda, lambda_net = 0.0;
+  element_lambda = malloc((cooling->N_Elements+2)*sizeof(double));
+
+  for (int j = 0; j < cooling->N_Elements+2; j++) element_lambda[j] = 0.0;
+  eagle_metal_cooling_rate(p, cooling, cosmo, internal_const, element_lambda);
+  for (int j = 0; j < cooling->N_Elements+2; j++) lambda_net += element_lambda[j];
+
+  return lambda_net;
+}
+
+__attribute__((always_inline)) INLINE static double eagle_print_metal_cooling_rate(const struct part* restrict p, const struct cooling_function_data* restrict cooling, const struct cosmology* restrict cosmo, const struct phys_const *internal_const) {
+  double *element_lambda, lambda_net = 0.0;
+  element_lambda = malloc((cooling->N_Elements+2)*sizeof(double));
+  
+  char output_filename[21];
+  FILE** output_file = malloc((cooling->N_Elements+2)*sizeof(FILE*));
+  for (int element = 0; element < cooling->N_Elements+2; element++){
+    sprintf(output_filename, "%s%d%s", "cooling_output_",element,".dat");
+    output_file[element] = fopen(output_filename, "a");
+    if (output_file == NULL)
+    {   
+        printf("Error opening file!\n");
+        exit(1);
+    }
+  }
+
+  for (int j = 0; j < cooling->N_Elements+2; j++) element_lambda[j] = 0.0;
+  eagle_metal_cooling_rate(p, cooling, cosmo, internal_const, element_lambda);
+  for (int j = 0; j < cooling->N_Elements+2; j++) {
+    lambda_net += element_lambda[j];
+    fprintf(output_file[j],"%.5e\n",element_lambda[j]);
+  }
+  
+  for (int i = 0; i < cooling->N_Elements+2; i++) fclose(output_file[i]);
 
-  return LambdaNet;
+  return lambda_net;
 }
 
 /**
@@ -651,6 +700,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
 
   XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
   HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] / (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
+  //printf("Eagle cooling.h density %.5e\n", hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY));
 
   /* convert Hydrogen mass fraction in Hydrogen number density */
   inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
@@ -845,7 +895,7 @@ static INLINE void cooling_init_backend(
   printf("Eagle cooling.h read table \n");
 
   cooling->internal_energy_scale = units_cgs_conversion_factor(us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
-  cooling->number_density_scale = units_cgs_conversion_factor(us,UNIT_CONV_DENSITY)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);// /(2.0*phys_const->const_proton_mass); // CHECK OUT NUMBER DENSITY SCALING!!!
+  cooling->number_density_scale = units_cgs_conversion_factor(us,UNIT_CONV_DENSITY)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
   cooling->power_scale = units_cgs_conversion_factor(us,UNIT_CONV_POWER)/units_cgs_conversion_factor(us,UNIT_CONV_MASS);
   cooling->temperature_scale = units_cgs_conversion_factor(us,UNIT_CONV_TEMPERATURE);
 

src/cooling/EAGLE/cooling_struct.h

@@ -71,6 +71,7 @@ struct cooling_tables {
  */
 struct eagle_cooling_table {
   struct cooling_tables_redshift_invariant photoionisation_cooling;
+  struct cooling_tables_redshift_invariant photodissociation_cooling;
   struct cooling_tables_redshift_invariant collisional_cooling;
   struct cooling_tables element_cooling;
 };

src/cooling/EAGLE/eagle_cool_tables.h

@@ -78,8 +78,6 @@ inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
   }
   fclose(infile);
 
-  //printf("eagle_cool_tables.h redshift max, min, N_Redshifts: %.5e, %.5e, %d\n", cooling->Redshifts[cooling->N_Redshifts-1], cooling->Redshifts[0], cooling->N_Redshifts);
-
 }
 
 inline void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
@@ -352,6 +350,111 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
 
   char fname[500], set_name[500];
 
+  int specs, i, j, table_index, cooling_index;
+
+  float *net_cooling_rate;
+  float *electron_abundance;
+  float *temperature;
+  float *he_net_cooling_rate;
+  float *he_electron_abundance;
+
+  net_cooling_rate = (float *)malloc(cooling->N_Temp*sizeof(float));
+  electron_abundance = (float *)malloc(cooling->N_Temp*sizeof(float));
+  temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+  he_net_cooling_rate = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+  he_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+
+  cooling_table.metal_heating = (float *)malloc(cooling->N_Elements*cooling->N_Temp*sizeof(float));
+  cooling_table.electron_abundance = (float *)malloc(cooling->N_Temp*sizeof(float));
+  cooling_table.temperature = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+  cooling_table.H_plus_He_heating = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+  cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*sizeof(float));
+
+  sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
+
+  file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
+
+  /* For normal elements */
+  for (specs = 0; specs < cooling->N_Elements; specs++) {
+    sprintf(set_name, "/%s/Net_Cooling", cooling->ElementNames[specs]);
+    dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+    status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                     net_cooling_rate);
+    status = H5Dclose(dataset);
+
+    for (j = 0; j < cooling->N_Temp; j++){
+        cooling_index = row_major_index_3d(0,specs,j,1,cooling->N_Elements,cooling->N_Temp); //Redshift invariant table!!!
+        cooling_table.metal_heating[cooling_index] = -net_cooling_rate[j];
+    }
+  }
+
+  /* Helium */
+  sprintf(set_name, "/Metal_free/Net_Cooling");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   he_net_cooling_rate);
+  status = H5Dclose(dataset);
+
+  sprintf(set_name, "/Metal_free/Temperature/Temperature");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   temperature);
+  status = H5Dclose(dataset);
+
+  sprintf(set_name, "/Metal_free/Electron_density_over_n_h");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   he_electron_abundance);
+  status = H5Dclose(dataset);
+
+  for (i = 0; i < cooling->N_He; i++){
+    for (j = 0; j < cooling->N_Temp; j++){
+      table_index = row_major_index_2d(i,j,cooling->N_He,cooling->N_Temp);
+      cooling_index = row_major_index_3d(0,i,j,1,cooling->N_He,cooling->N_Temp); //Redshift invariant table!!!
+      cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
+      cooling_table.H_plus_He_electron_abundance[cooling_index] =
+          he_electron_abundance[table_index];
+      cooling_table.temperature[cooling_index] = log10(temperature[table_index]);
+    }
+  }
+
+  sprintf(set_name, "/Solar/Electron_density_over_n_h");
+  dataset = H5Dopen(file_id, set_name, H5P_DEFAULT);
+  status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
+                   electron_abundance);
+  status = H5Dclose(dataset);
+
+  for (i = 0; i < cooling->N_Temp; i++){
+      cooling_table.electron_abundance[i] = electron_abundance[i];
+  }
+
+  status = H5Fclose(file_id);
+
+  free(net_cooling_rate);
+  free(electron_abundance);
+  free(temperature);
+  free(he_net_cooling_rate);
+  free(he_electron_abundance);
+
+  printf("eagle_cool_tables.h done reading in collisional table\n");
+  return cooling_table;
+}
+
+/*
+ * ----------------------------------------------------------------------
+ * Get the cooling table for photodissociation 
+ * ----------------------------------------------------------------------
+ */
+
+inline struct cooling_tables_redshift_invariant get_photodis_table(char *cooling_table_path, const struct cooling_function_data* restrict cooling) {
+
+  struct cooling_tables_redshift_invariant cooling_table;
+  hid_t file_id, dataset;
+
+  herr_t status;
+
+  char fname[500], set_name[500];
+
   int specs, i, j, k, table_index, cooling_index;
 
   float *net_cooling_rate;
@@ -373,7 +476,6 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
   cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
 
   sprintf(fname, "%sz_photodis.hdf5", cooling_table_path);
-  //sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
 
   file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
 
@@ -388,7 +490,7 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
     for (j = 0; j < cooling->N_Temp; j++){
       for (k = 0; k < cooling->N_nH; k++){
         table_index = row_major_index_2d(j,k,cooling->N_Temp,cooling->N_nH);
-        cooling_index = row_major_index_4d(0,specs,k,j,1,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
+        cooling_index = row_major_index_3d(specs,j,k,cooling->N_Elements,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
         cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
       }
     }
@@ -417,7 +519,7 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
     for (j = 0; j < cooling->N_Temp; j++){
       for (k = 0; k < cooling->N_nH; k++) {
         table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
-        cooling_index = row_major_index_4d(0,i,k,j,1,cooling->N_He,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
+        cooling_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
         cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
         cooling_table.H_plus_He_electron_abundance[cooling_index] =
             he_electron_abundance[table_index];
@@ -435,26 +537,20 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
   for (i = 0; i < cooling->N_Temp; i++){
     for (j = 0; j < cooling->N_nH; j++){
       table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
-      cooling_index = row_major_index_3d(0,j,i,1,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
+      cooling_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH); //Redshift invariant table!!!
       cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
     }
   }
 
   status = H5Fclose(file_id);
 
-  //cooling_table.metal_heating = cooling_MetalsNetHeating;
-  //cooling_table.H_plus_He_heating = cooling_HplusHeNetHeating;
-  //cooling_table.H_plus_He_electron_abundance = cooling_HplusHeElectronAbundance;
-  //cooling_table.temperature = cooling_ThermalToTemp;
-  //cooling_table.electron_abundance = cooling_SolarElectronAbundance;
-
   free(net_cooling_rate);
   free(electron_abundance);
   free(temperature);
   free(he_net_cooling_rate);
   free(he_electron_abundance);
 
-  printf("eagle_cool_tables.h done reading in collisional table\n");
+  printf("eagle_cool_tables.h done reading in photodissociation table\n");
   return cooling_table;
 }
 
@@ -513,7 +609,8 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
       for (i = 0; i < cooling->N_nH; i++){
         for (j = 0; j < cooling->N_Temp; j++){
           table_index = row_major_index_2d(j,i,cooling->N_Temp,cooling->N_nH);
-          cooling_index = row_major_index_4d(z_index,specs,i,j,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); 
+          //cooling_index = row_major_index_4d(z_index,specs,i,j,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); 
+          cooling_index = row_major_index_3d(specs,i,j,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); 
           cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
         }
       }
@@ -542,7 +639,8 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
       for (j = 0; j < cooling->N_Temp; j++){
         for (k = 0; k < cooling->N_nH; k++) {
           table_index = row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
-          cooling_index = row_major_index_4d(z_index,i,k,j,cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp); 
+          //cooling_index = row_major_index_4d(z_index,i,k,j,cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp); 
+          cooling_index = row_major_index_3d(i,k,j,cooling->N_He,cooling->N_nH,cooling->N_Temp); 
           cooling_table.H_plus_He_heating[cooling_index] = -he_net_cooling_rate[table_index];
           cooling_table.H_plus_He_electron_abundance[cooling_index] =
               he_electron_abundance[table_index];
@@ -560,7 +658,8 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
     for (i = 0; i < cooling->N_Temp; i++){
       for (j = 0; j < cooling->N_nH; j++){
         table_index = row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
-        cooling_index = row_major_index_3d(z_index,j,i,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); 
+        //cooling_index = row_major_index_3d(z_index,j,i,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); 
+        cooling_index = row_major_index_2d(j,i,cooling->N_nH,cooling->N_Temp); 
         cooling_table.electron_abundance[cooling_index] = electron_abundance[table_index];
       }
     }
@@ -568,12 +667,6 @@ inline struct cooling_tables get_cooling_table(char *cooling_table_path, const s
     status = H5Fclose(file_id);
   }
 
-  //cooling_table.metal_heating = cooling_MetalsNetHeating;
-  //cooling_table.H_plus_He_heating = cooling_HplusHeNetHeating;
-  //cooling_table.H_plus_He_electron_abundance = cooling_HplusHeElectronAbundance;
-  //cooling_table.temperature = cooling_ThermalToTemp;
-  //cooling_table.electron_abundance = cooling_SolarElectronAbundance;
-  
   free(net_cooling_rate);
   free(electron_abundance);
   free(temperature);
@@ -590,6 +683,7 @@ inline struct eagle_cooling_table eagle_readtable(char *cooling_table_path, cons
   struct eagle_cooling_table table;
 
   table.photoionisation_cooling = get_no_compt_table(cooling_table_path, cooling);
+  table.photodissociation_cooling = get_photodis_table(cooling_table_path, cooling);
   table.collisional_cooling = get_collisional_table(cooling_table_path, cooling);
   table.element_cooling = get_cooling_table(cooling_table_path, cooling);
 
@@ -662,107 +756,4 @@ inline void MakeNamePointers(struct cooling_function_data* cooling) {
   }
 }
 
-
-//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
+elements = 11
 temperature = []
 cooling_rate = [[] for i in range(elements+1)]
 length = 0
@@ -21,18 +21,20 @@ for elem in range(elements):
 
 
 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')
+p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
+p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
+p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
+p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
+p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p10, = plt.loglog(temperature, cooling_rate[11], 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.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
 plt.show()

tests/testCooling.c

@@ -95,7 +95,7 @@ int main() {
     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));
 
-    cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const);
+    cooling_du_dt = eagle_print_metal_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",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);