diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 4bc6ff79736c0c1f5ebacde6b642ff13a9fe84ed..b468cb059e559392d57a7673d069d1ac3781f607 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/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)
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 8d73e215044a034148981d4877e036e01a15d1e4..23de7674c8fddfc44e6f5201b4f35f3dc5dc808e 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/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
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 83373ca99e656d5adc72d64977042414635661aa..232bd54650330dafc62d37bbd97552ed319ad732 100755
--- a/examples/CoolingBox/run.sh
+++ b/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
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index b542f9e99ec365107a5c4f265e91bb4f13b1e844..e2c5b7e811d39c6e3bc4ba2cb8c8ac5a62108b16 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/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);
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index df63c5a2a5cc5ecaf946e2f4f93ff99f050edbad..b9f587108eb06d135ecfa72fa3fdf81ed2b3f1e5 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/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;
};
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index 1c9f068b0d2dcc8efc82a9573e3df62a8eef6333..1b0b33aa6f70494a4277297e37a47b74c99b3fd3 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/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
diff --git a/tests/plots.py b/tests/plots.py
index 844dea6d1745b99650b96d7a3a8d4de1a4936665..56d90d1dc79fed9c6ca9e24b64697c45a30f9baa 100644
--- a/tests/plots.py
+++ b/tests/plots.py
@@ -1,6 +1,6 @@
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()
diff --git a/tests/testCooling.c b/tests/testCooling.c
index be35f7969b907358665a3d3ca23691bed9e617c5..62b05dc895c365d3f4bb23497bc8461cfe7bd503 100644
--- a/tests/testCooling.c
+++ b/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);