diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
new file mode 100644
index 0000000000000000000000000000000000000000..4bc6ff79736c0c1f5ebacde6b642ff13a9fe84ed
--- /dev/null
+++ b/examples/CoolingBox/analytical_test.py
@@ -0,0 +1,149 @@
+import matplotlib
+matplotlib.use("Agg")
+from pylab import *
+import h5py
+
+# Plot parameters
+params = {'axes.labelsize': 10,
+'axes.titlesize': 10,
+'font.size': 12,
+'legend.fontsize': 12,
+'xtick.labelsize': 10,
+'ytick.labelsize': 10,
+'text.usetex': True,
+ 'figure.figsize' : (3.15,3.15),
+'figure.subplot.left' : 0.2,
+'figure.subplot.right' : 0.99,
+'figure.subplot.bottom' : 0.2,
+'figure.subplot.top' : 0.9,
+'figure.subplot.wspace' : 0.15,
+'figure.subplot.hspace' : 0.12,
+'lines.markersize' : 6,
+'lines.linewidth' : 3.,
+'text.latex.unicode': True
+}
+rcParams.update(params)
+rc('font',**{'family':'sans-serif','sans-serif':['Times']})
+
+
+import numpy as np
+import h5py as h5
+import sys
+
+def interpol_lambda(temp_list,cooling_rate,temp):
+ #print(temp,temp_list[0],temp_list[len(temp_list)-1])
+ if temp < temp_list[0]:
+ return[cooling_rate[0]]
+ if temp > temp_list[len(temp_list)-1]:
+ return[cooling_rate[len(cooling_rate)-1]]
+ j = 0
+ while temp_list[j+1] < temp:
+ j += 1
+ cooling = cooling_rate[j]*(temp_list[j+1]-temp)/(temp_list[j+1]-temp_list[j]) + cooling_rate[j+1]*(temp-temp_list[j])/(temp_list[j+1]-temp_list[j])
+ return cooling
+
+def convert_u_to_temp_sol(u,rho):
+ k_b = 1.38E-16 #boltzmann
+ m_p = 1.67e-24 #proton mass
+ gamma = 5.0/3.0
+ pressure = u*rho*(gamma - 1.0)
+ temp = pressure/(k_b*rho/(0.59*m_p))
+ return temp
+
+# File containing the total energy
+stats_filename = "./energy.txt"
+
+# First snapshot
+snap_filename = "coolingBox_0000.hdf5"
+
+# Some constants in cgs units
+k_b = 1.38E-16 #boltzmann
+m_p = 1.67e-24 #proton mass
+
+# Initial conditions set in makeIC.py
+T_init = 1.0e7
+
+# Read the initial state of the gas
+f = h5.File(snap_filename,'r')
+rho = np.mean(f["/PartType0/Density"])
+pressure = np.mean(f["/PartType0/Pressure"])
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+unit_vel = unit_length/unit_time
+#hubble_param = 0.71
+hubble_param = 1.0
+
+unit_length = unit_length/hubble_param
+unit_mass = unit_mass/hubble_param
+
+yHe = 0.28
+temp_0 = 1.0e7
+
+rho = rho*unit_mass/(unit_length**3)
+
+# Read snapshots
+nsnap = 40
+npart = 32768
+u_snapshots_cgs = zeros(nsnap)
+u_part_snapshots_cgs = zeros((nsnap,npart))
+t_snapshots_cgs = zeros(nsnap)
+for i in range(nsnap):
+ snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
+ u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*unit_length**2/(unit_time**2)
+ t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
+
+# calculate analytic solution. Since cooling rate is constant,
+# temperature decrease in linear in time.
+# read Lambda and temperatures from table
+temperature = []
+cooling_rate = []
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+ data = line.split()
+ temperature.append(float(data[0]))
+ cooling_rate.append(-float(data[1]))
+
+tfinal = 3.3*t_snapshots_cgs[nsnap-1]
+nt = 1e4
+dt = tfinal/nt
+
+t_sol = np.zeros(int(nt))
+temp_sol = np.zeros(int(nt))
+lambda_sol = np.zeros(int(nt))
+u = np.mean(u_part_snapshots_cgs[0,:])
+temp_sol[0] = convert_u_to_temp_sol(u,rho)
+#print(u,temp_sol[0])
+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)
+ 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)
+ temp_sol[j+1] = convert_u_to_temp_sol(u_next,rho)
+ lambda_sol[j] = Lambda_net
+ u = u_next
+
+mean_temp = np.zeros(nsnap)
+for j in range(nsnap):
+ mean_temp[j] = convert_u_to_temp_sol(np.mean(u_part_snapshots_cgs[j,:]),rho)
+p = figure()
+p1, = plot(t_sol,temp_sol,linewidth = 0.5,color = 'k',label = 'analytical')
+p2, = plot(t_snapshots_cgs,mean_temp,linewidth = 0.5,color = 'r',label = 'swift mean')
+l = legend(handles = [p1,p2])
+xlabel("${\\rm{Time~[s]}}$", labelpad=0)
+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}}$")
+
+savefig("cooling.png", dpi=200)
diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index b4af701343ae309bee2d9858f8620e1d4d79a44d..720179b767dcf65d60e57ca88d418b278d889d21 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/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.25 # The end time of the simulation (in internal units).
+ time_end: 0.02 # 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: 1e-2 # Time difference between consecutive outputs (in internal units)
+ delta_time: 5e-4 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index f863e174b1fcd404ae178fe324c7a165598b4af0..8d73e215044a034148981d4877e036e01a15d1e4 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.5e6 # Pressure in code units (at 10^5K)
+P = 4.5e8 # Pressure in code units (at 10^7K)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 30b2177a6e8bb95a20146397f8b6a5021161b27f..83373ca99e656d5adc72d64977042414635661aa 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -21,7 +21,7 @@ then
fi
# Run SWIFT
-../swift -s -C -t 1 coolingBox.yml
+../swift -s -C -t 4 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 57e2d313edea5b76b013b0c55e1ea2c2bf511658..b542f9e99ec365107a5c4f265e91bb4f13b1e844 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -41,7 +41,6 @@
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
-// #include "physical_constants_cgs.h"
#include "units.h"
#include "eagle_cool_tables.h"
@@ -111,6 +110,8 @@ __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;
}
}
@@ -129,7 +130,7 @@ __attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int
*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);
+ //printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
*dx = ((float)x - table[*i]) * dxm1;
}
}
@@ -184,6 +185,8 @@ __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;
}
@@ -318,7 +321,8 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int
inline int set_cooling_SolarAbundances(const float *element_abundance,
double *cooling_element_abundance,
- const struct cooling_function_data* restrict cooling) {
+ const struct cooling_function_data* restrict cooling,
+ const struct part* restrict p) {
int i, index;
int Silicon_SPH_Index = -1;
int Calcium_SPH_Index = -1;
@@ -372,30 +376,56 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
if (strcmp(chemistry_get_element_name((enum chemistry_element) i), "Silicon") == 0) sili_index = i;
}
+ // Eagle way of identifying and assigning element abundance with strange workaround for calcium and sulphur
+ //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[cooling->ElementNamePointers[i]] *
+ // cooling_ElementAbundance_SOLARM1[i];
+ // //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
+ // }
+ //}
+
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 (i == Calcium_CoolHeat_Index && Calcium_SPH_Index != -1)
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];
+ element_abundance[Silicon_SPH_Index] * cooling->calcium_over_silicon_ratio *
+ cooling_ElementAbundance_SOLARM1[Calcium_CoolHeat_Index];
}
- else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
+ 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];
+ element_abundance[Silicon_SPH_Index] * cooling->sulphur_over_silicon_ratio *
+ cooling_ElementAbundance_SOLARM1[Sulphur_CoolHeat_Index];
}
else{
cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
cooling_ElementAbundance_SOLARM1[i];
- //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
}
+ //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i], element_abundance[i]);
}
return 0;
@@ -475,10 +505,13 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons
get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
+ //logT = interpol_2d(cooling->table.collisional_cooling.temperature, u_i, He_i, d_u, d_He, cooling->N_Temp, cooling->N_He);
+
logT = interpol_4d(cooling->table.element_cooling.temperature, 0, He_i, n_h_i,
u_i, dz, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
+
T = pow(10.0, logT);
if (u_i == 0 && d_u == 0) T *= u / pow(10.0, cooling->Therm[0]);
@@ -530,10 +563,10 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* ------------------ */
LambdaNet =
- interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, 0, He_i, n_h_i,
+ 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);
h_plus_he_electron_abundance =
- interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
+ 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);
/* ------------------ */
@@ -558,15 +591,15 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* for each element, find the abundance and multiply it
by the interpolated heating-cooling */
- set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling);
+ set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling, p);
solar_electron_abundance =
- interpol_3d(cooling->table.element_cooling.electron_abundance, 0, n_h_i, temp_i, dz,
+ 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 */
double element_cooling_lambda;
for (i = 0; i < cooling->N_Elements; i++){
- element_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_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) *
(h_plus_he_electron_abundance / solar_electron_abundance) *
cooling_solar_abundances[i];
@@ -610,6 +643,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
int i;
static double zold = 100, LambdaCumul = 0;
+ dt *= units_cgs_conversion_factor(us,UNIT_CONV_TIME);
u = u_old;
u_lower = u;
@@ -620,6 +654,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* 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;
+ //inn_h = hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY) * XH / eagle_proton_mass_cgs;
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
@@ -707,16 +742,12 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
}
float cooling_du_dt = 0.0;
- if (dt > 0) cooling_du_dt = (u - u_old)/dt/cooling->power_scale;
- const float hydro_du_dt = hydro_get_internal_energy_dt(p);
-
- /* Integrate cooling equation to enforce energy floor */
- if (u_old + hydro_du_dt * dt < cooling->min_energy) {
- cooling_du_dt = 0.f;
- } else if (u_old + (hydro_du_dt + cooling_du_dt) * dt < cooling->min_energy) {
- cooling_du_dt = (u_old + dt * hydro_du_dt - cooling->min_energy) / dt;
+ if (dt > 0){
+ cooling_du_dt = (u - u_old)/dt/cooling->power_scale;
}
+ const float hydro_du_dt = hydro_get_internal_energy_dt(p);
+
/* Update the internal energy time derivative */
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
@@ -804,6 +835,8 @@ static INLINE void cooling_init_backend(
parser_get_param_string(parameter_file, "EagleCooling:filename",cooling->cooling_table_path);
cooling->reionisation_redshift = parser_get_param_float(parameter_file, "EagleCooling:reionisation_redshift");
+ cooling->calcium_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:CalciumOverSilicon");
+ cooling->sulphur_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:SulphurOverSilicon");
GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname,cooling);
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index e9f488549bb0c3c45c5775fdebbc6f8ad9c3f506..df63c5a2a5cc5ecaf946e2f4f93ff99f050edbad 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -113,6 +113,8 @@ struct cooling_function_data {
double power_scale;
char cooling_table_path[500];
float reionisation_redshift;
+ float calcium_over_silicon_ratio;
+ float sulphur_over_silicon_ratio;
};
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index b25496f379f9a1b95740ab46fd0a5f9fa203985e..1c9f068b0d2dcc8efc82a9573e3df62a8eef6333 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -261,24 +261,15 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
net_cooling_rate);
status = H5Dclose(dataset);
- //for(i = 0; i < cooling->N_Temp*cooling->N_nH; i++) printf("eagle_cool_tables.h i, net_cooling_rate %d, %.5e\n",i,net_cooling_rate[i]);
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_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
- //printf("eagle_cool_tables.h j,k,j*cooling->N_nH+k,table_index,cooling_index,net cooling rate, cooling_table value %d, %d, %d, %d, %d %.5e, %.5e\n",j,k,j*cooling->N_nH+k,table_index,cooling_index,-net_cooling_rate[table_index],cooling_table.metal_heating[cooling_index]);
}
}
}
- //for (i = 0; i < cooling->N_nH; i++){
- // table_index = row_major_index_2d(0,i,cooling->N_Temp,cooling->N_nH);
- // cooling_index = row_major_index_4d(0,0,i,0,1,cooling->N_Elements,cooling->N_nH,cooling->N_Temp);
- // printf("eagle_cool_tables.h i, cooling_index, table_index, cooling_table.metal_heating, net heating = %d %d %d %.5e %.5e\n",i,cooling_index,table_index,cooling_table.metal_heating[cooling_index],net_cooling_rate[table_index]);
- //}
- //for (i = 0; i < cooling->N_Elements*cooling->N_Temp*cooling->N_nH; i++) printf("eagle cooling.h i, cooling_table %d, %.5e\n",i,cooling_table.metal_heating[i]);
- //error("stop in eagle_cool_tables.h");
/* Helium */
@@ -382,6 +373,7 @@ 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);
diff --git a/src/cooling/EAGLE/eagle_init_cool.h b/src/cooling/EAGLE/eagle_init_cool.h
deleted file mode 100644
index 7081d45b1bac90bbbf6ebafd9c298da30a383560..0000000000000000000000000000000000000000
--- a/src/cooling/EAGLE/eagle_init_cool.h
+++ /dev/null
@@ -1,899 +0,0 @@
-#ifndef SWIFT_EAGLE_INIT_COOL_H
-#define SWIFT_EAGLE_INIT_COOL_H
-
-
-/*
- * ----------------------------------------------------------------------
- * Prototype for built-in testing routines
- * ----------------------------------------------------------------------
- */
-
-void DriverEagleTestCool();
-
-void EagleTestCool(float z);
-
-void DriverEagleTestCoolThermalEvolution();
-
-/*
- * ----------------------------------------------------------------------
- * This routine gets the redshifts from a text file
- * ----------------------------------------------------------------------
- */
-
-void GetCoolingRedshifts() {
- FILE *infile;
-
- int i = 0;
-
- char buffer[500], redfilename[500];
-
- sprintf(redfilename, "%s/redshifts.dat", All.CoolTablePath);
- infile = fopen(redfilename, "r");
- if (infile == NULL) puts("GetCoolingRedshifts can't open a file");
-
- if (fscanf(infile, "%s", buffer) != EOF) {
- cooling_N_Redshifts = atoi(buffer);
- cooling_Redshifts =
- (float *)mymalloc("cNz", cooling_N_Redshifts * sizeof(float));
-
- while (fscanf(infile, "%s", buffer) != EOF) {
- cooling_Redshifts[i] = atof(buffer);
- i += 1;
- }
- }
- fclose(infile);
-
-#ifdef EAGLE_VERBOSE
- for (i = 0; i < cooling_N_Redshifts; i++)
- printf("Cooling table %2d at redshift %1.3f\n", i + 1,
- cooling_Redshifts[i]);
- printf("\n");
-#endif
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine reads in the header of the cooling table files
- * ----------------------------------------------------------------------
- */
-
-void ReadCoolingHeader(char *fname) {
- int i;
-
- hid_t tempfile_id, dataset, datatype;
-
- herr_t status;
-
- /* fill the constants */
- tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
- if (tempfile_id < 0) {
- printf("[ReadCoolingHeader()]: unable to open file %s\n", fname);
- endrun(101);
- }
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_density_bins");
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling_N_nH);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins");
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling_N_Temp);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions");
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling_N_He);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals");
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling_N_Elements);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances");
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling_N_SolarAbundances);
- status = H5Dclose(dataset);
-
- /* allocate arrays for cooling table header */
- allocate_header_arrays();
-
- /* fill the arrays */
- dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins");
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling_Temp);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Solar/Hydrogen_density_bins");
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling_nH);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins");
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling_Therm);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins");
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling_HeFrac);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions");
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling_SolarAbundances);
- status = H5Dclose(dataset);
-
- char element_names[cooling_N_Elements][EL_NAME_LENGTH];
- hsize_t string_length = EL_NAME_LENGTH;
-
- /* names of chemical elements stored in table */
- datatype = H5Tcopy(H5T_C_S1);
- status = H5Tset_size(datatype, string_length);
- dataset = H5Dopen(tempfile_id, "/Header/Metal_names");
- status =
- H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
- status = H5Dclose(dataset);
-
- for (i = 0; i < cooling_N_Elements; i++)
- cooling_ElementNames[i] = mystrdup(element_names[i]);
-
- char solar_abund_names[cooling_N_SolarAbundances][EL_NAME_LENGTH];
-
- /* assumed solar abundances used in constructing the tables, and corresponding
- * names */
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names");
- status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- solar_abund_names);
- status = H5Dclose(dataset);
- H5Tclose(datatype);
-
- for (i = 0; i < cooling_N_SolarAbundances; i++)
- cooling_SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
-
- status = H5Fclose(tempfile_id);
-
- /* Convert to temperature, density and internal energy arrays to log10 */
- for (i = 0; i < cooling_N_Temp; i++) {
- cooling_Temp[i] = log10(cooling_Temp[i]);
- cooling_Therm[i] = log10(cooling_Therm[i]);
- }
-
- for (i = 0; i < cooling_N_nH; i++) cooling_nH[i] = log10(cooling_nH[i]);
-
- printf("Done with cooling table header.\n");
- fflush(stdout);
-
-#ifdef EAGLE_VERBOSE
- printf(
- "\n\nTemperature grid (LOG10) for energy to temperature conversion\n\n");
- for (i = 0; i < cooling_N_Temp; i++)
- printf("Temperature[%3d] = %f\n", i, cooling_Temp[i]);
-
- printf("\n\nEnergy grid (LOG10) for energy to temperature conversion\n\n");
- for (i = 0; i < cooling_N_Temp; i++)
- printf("Energy[%3d] = %f\n", i, cooling_Therm[i]);
-
- printf("\n\nDensity grid (LOG10) for energy to temperature conversion\n\n");
- for (i = 0; i < cooling_N_nH; i++)
- printf("Density[%2d] = %f\n", i, cooling_nH[i]);
-#endif
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine broadcasts the header infos
- * ----------------------------------------------------------------------
- */
-
-void BroadCastHeader() {
- int position, bufsize, size, i;
-
- char *buffer;
-
- /* get size of the buffer */
- MPI_Pack_size(6, MPI_INT, MPI_COMM_WORLD, &size);
- bufsize = size;
-
- /* allocate memory for the buffer */
- buffer = (char *)mymalloc("buffer", bufsize);
-
- if (ThisTask == 0) {
- position = 0;
-
- /* pack the array dimensions */
- MPI_Pack(&cooling_N_nH, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(&cooling_N_Temp, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(&cooling_N_He, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(&cooling_N_Elements, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(&cooling_N_SolarAbundances, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(&cooling_N_Redshifts, 1, MPI_INT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- }
-
- MPI_Bcast(buffer, bufsize, MPI_PACKED, 0, MPI_COMM_WORLD);
-
- if (ThisTask != 0) {
- position = 0;
-
- /* unpack the array dimensions */
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_nH, 1, MPI_INT,
- MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_Temp, 1, MPI_INT,
- MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_He, 1, MPI_INT,
- MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_Elements, 1, MPI_INT,
- MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_SolarAbundances, 1,
- MPI_INT, MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, &cooling_N_Redshifts, 1, MPI_INT,
- MPI_COMM_WORLD);
- }
-
- /* free allocated memory */
- myfree(buffer);
-
- if (ThisTask != 0) {
- /* allocate header arrays */
- allocate_header_arrays();
-
- /* allocate redshift array */
- cooling_Redshifts =
- (float *)mymalloc("cNz", cooling_N_Redshifts * sizeof(float));
- }
-
- bufsize = 0;
- /* get size of the buffer */
- MPI_Pack_size(cooling_N_Temp, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_nH, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_Temp, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_He, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_SolarAbundances, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_Redshifts, MPI_FLOAT, MPI_COMM_WORLD, &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_Elements * EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD,
- &size);
- bufsize += size;
- MPI_Pack_size(cooling_N_SolarAbundances * EL_NAME_LENGTH, MPI_CHAR,
- MPI_COMM_WORLD, &size);
- bufsize += size;
-
- /* allocate memory for the buffer */
- buffer = (char *)mymalloc("buffer", bufsize);
-
- if (ThisTask == 0) {
- position = 0;
-
- MPI_Pack(cooling_Temp, cooling_N_Temp, MPI_FLOAT, buffer, bufsize,
- &position, MPI_COMM_WORLD);
- MPI_Pack(cooling_nH, cooling_N_nH, MPI_FLOAT, buffer, bufsize, &position,
- MPI_COMM_WORLD);
- MPI_Pack(cooling_Therm, cooling_N_Temp, MPI_FLOAT, buffer, bufsize,
- &position, MPI_COMM_WORLD);
- MPI_Pack(cooling_HeFrac, cooling_N_He, MPI_FLOAT, buffer, bufsize,
- &position, MPI_COMM_WORLD);
- MPI_Pack(cooling_SolarAbundances, cooling_N_SolarAbundances, MPI_FLOAT,
- buffer, bufsize, &position, MPI_COMM_WORLD);
- MPI_Pack(cooling_Redshifts, cooling_N_Redshifts, MPI_FLOAT, buffer, bufsize,
- &position, MPI_COMM_WORLD);
-
- for (i = 0; i < cooling_N_Elements; i++)
- MPI_Pack(cooling_ElementNames[i], EL_NAME_LENGTH, MPI_CHAR, buffer,
- bufsize, &position, MPI_COMM_WORLD);
- for (i = 0; i < cooling_N_SolarAbundances; i++)
- MPI_Pack(cooling_SolarAbundanceNames[i], EL_NAME_LENGTH, MPI_CHAR, buffer,
- bufsize, &position, MPI_COMM_WORLD);
- }
-
- MPI_Bcast(buffer, bufsize, MPI_PACKED, 0, MPI_COMM_WORLD);
-
- if (ThisTask != 0) {
- position = 0;
-
- MPI_Unpack(buffer, bufsize, &position, cooling_Temp, cooling_N_Temp,
- MPI_FLOAT, MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, cooling_nH, cooling_N_nH, MPI_FLOAT,
- MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, cooling_Therm, cooling_N_Temp,
- MPI_FLOAT, MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, cooling_HeFrac, cooling_N_He,
- MPI_FLOAT, MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, cooling_SolarAbundances,
- cooling_N_SolarAbundances, MPI_FLOAT, MPI_COMM_WORLD);
- MPI_Unpack(buffer, bufsize, &position, cooling_Redshifts,
- cooling_N_Redshifts, MPI_FLOAT, MPI_COMM_WORLD);
-
- for (i = 0; i < cooling_N_Elements; i++)
- MPI_Unpack(buffer, bufsize, &position, cooling_ElementNames[i],
- EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD);
- for (i = 0; i < cooling_N_SolarAbundances; i++)
- MPI_Unpack(buffer, bufsize, &position, cooling_SolarAbundanceNames[i],
- EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD);
- }
-
- /* free allocated memory */
- myfree(buffer);
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine sets the value of solar metallicity using solar element
- * abundances stored in the cooling tables.
- * ----------------------------------------------------------------------
- */
-
-void set_solar_metallicity(void) {
- int i;
-
- if (ThisTask == 0) printf("\n\nSolar abundances:\n\n");
-
- /* set solar metallicity value */
- for (i = 0, SolarMetallicity = 1; i < cooling_N_SolarAbundances; i++) {
- if (ThisTask == 0)
- printf("Solar abundance of %s = %g\n", cooling_SolarAbundanceNames[i],
- cooling_SolarAbundances[i]);
- if (strcmp(cooling_SolarAbundanceNames[i], "Hydrogen") == 0)
- SolarMetallicity -= cooling_SolarAbundances[i];
-
- if (strcmp(cooling_SolarAbundanceNames[i], "Helium") == 0)
- SolarMetallicity -= cooling_SolarAbundances[i];
- }
- if (ThisTask == 0) printf("\nSolarMetallicity = %g\n", SolarMetallicity);
-}
-
-/*
- * ----------------------------------------------------------------------
- * Make pointers to element names
- * ----------------------------------------------------------------------
- */
-
-void MakeNamePointers() {
- int i, j, sili_index = 0;
-
- ElementNamePointers =
- (int *)mymalloc("coolNele", cooling_N_Elements * sizeof(int));
- SolarAbundanceNamePointers =
- (int *)mymalloc("coolNele2", cooling_N_Elements * sizeof(int));
-
- for (i = 0; i < EAGLE_NELEMENTS; i++) {
- if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
- }
-
- for (i = 0; i < cooling_N_Elements; i++) {
- SolarAbundanceNamePointers[i] = -999;
- ElementNamePointers[i] = -999;
-
- for (j = 0; j < cooling_N_SolarAbundances; j++) {
- if (strcmp(cooling_ElementNames[i], cooling_SolarAbundanceNames[j]) == 0)
- SolarAbundanceNamePointers[i] = j;
- }
-
- if (strcmp(cooling_ElementNames[i], "Sulphur") == 0 ||
- strcmp(cooling_ElementNames[i], "Calcium") ==
- 0) /* These elements are tracked! */
- ElementNamePointers[i] = -1 * sili_index;
- else {
- for (j = 0; j < EAGLE_NELEMENTS; j++) {
- if (strcmp(cooling_ElementNames[i], ElementNames[j]) == 0)
- ElementNamePointers[i] = j;
- }
- }
- }
-
-#ifdef EAGLE_VERBOSE
- if (ThisTask == 0) {
- for (i = 0; i < cooling_N_Elements; i++)
- printf("cooling_ElementNames[%d] is cooling_SolarAbundanceNames %s\n", i,
- cooling_SolarAbundanceNames[SolarAbundanceNamePointers[i]]);
- printf("\n");
-
- for (i = 0; i < cooling_N_Elements; i++)
- printf("cooling_ElementNames[%d] is ElementNames %s\n", i,
- ElementNames[ElementNamePointers[i]]);
- printf("\n");
- }
-#endif
-}
-
-void eagle_set_abundances_to_solar() {
- int i, element, index;
- if (ThisTask == 0)
- printf(" setting abundances of stars and particles to solar \n");
-
- for (element = 0; element < cooling_N_Elements; element++) {
- index = get_element_index(cooling_SolarAbundanceNames,
- cooling_N_SolarAbundances, ElementNames[element]);
- if (index >= 0) {
- for (i = 0; i < NumPart; i++)
- if (P[i].Type == 0) {
- SphP[i].Metals[element] = cooling_SolarAbundances[element];
- SphP[i].Metallicity = SolarMetallicity;
- } else if (P[i].Type == 4)
- StarP[P[i].pt.StarID].Metals[element] =
- cooling_SolarAbundances[element];
- }
- }
-}
-/*
- * ----------------------------------------------------------------------
- * This routine initialize cooling at the beginning of the simulation
- * ----------------------------------------------------------------------
- */
-void eagle_init_cool() {
- int i, j, k;
-
- double Z_Solar_Calcium = 0, Z_Solar_Sulphur = 0, Z_Solar_Silicon = 0;
-
- char fname[500];
-
- /* use a general file to get the header */
- if (ThisTask == 0) {
- GetCoolingRedshifts();
- sprintf(fname, "%sz_0.000.hdf5", All.CoolTablePath);
- ReadCoolingHeader(fname);
- }
-
- //BroadCastHeader();
-
- set_solar_metallicity();
-
- if (element_present("Silicon") == 0) {
- for (i = 0; i < cooling_N_SolarAbundances; i++) {
- if (strcmp(cooling_SolarAbundanceNames[i], "Calcium") == 0)
- Z_Solar_Calcium = cooling_SolarAbundances[i];
-
- if (strcmp(cooling_SolarAbundanceNames[i], "Sulphur") == 0)
- Z_Solar_Sulphur = cooling_SolarAbundances[i];
-
- if (strcmp(cooling_SolarAbundanceNames[i], "Silicon") == 0)
- Z_Solar_Silicon = cooling_SolarAbundances[i];
- }
- } else {
- printf("Silicon is not present: Calcium and Sulphur will not be tracked!");
- }
-
- //cooling_ThermalToTemp =
- // (float ****)mymalloc("coolTtoT", 2 * sizeof(float ***));
-
- //for (i = 0; i < 2; i++) {
- // cooling_ThermalToTemp[i] =
- // (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
-
- // for (j = 0; j < cooling_N_He; j++) {
- // cooling_ThermalToTemp[i][j] =
- // (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
-
- // for (k = 0; k < cooling_N_nH; k++)
- // cooling_ThermalToTemp[i][j][k] =
- // (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
- // }
- //}
-
- float cooling_ThermalToTemp[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
-
- //cooling_SolarElectronAbundance =
- // (float ***)mymalloc("coolSe", 2 * sizeof(float **));
-
- //for (i = 0; i < 2; i++) {
- // cooling_SolarElectronAbundance[i] =
- // (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
-
- // for (j = 0; j < cooling_N_nH; j++)
- // cooling_SolarElectronAbundance[i][j] =
- // (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
- //}
-
- float cooling_SolarElectronAbundance[2][cooling_N_nH][cooling_N_Temp];
-
- /* Allocate for metal cooling */
- //cooling_MetalsNetHeating =
- // (float ****)mymalloc("CoolHeat", 2 * sizeof(float ***));
-
- //for (i = 0; i < 2; i++) {
- // cooling_MetalsNetHeating[i] =
- // (float ***)mymalloc("coolHeat", cooling_N_Elements * sizeof(float **));
-
- // for (j = 0; j < cooling_N_Elements; j++) {
- // cooling_MetalsNetHeating[i][j] =
- // (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
-
- // for (k = 0; k < cooling_N_nH; k++)
- // cooling_MetalsNetHeating[i][j][k] =
- // (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
- // }
- //}
-
- float cooling_MetalsNetHeating[2][cooling_N_Elements][cooling_N_nH][cooling_N_Temp];
-
- /* Allocate for H and He cooling */
- //cooling_HplusHeNetHeating =
- // (float ****)mymalloc("CoolHHeHeat", 2 * sizeof(float ***));
-
- //for (i = 0; i < 2; i++) {
- // cooling_HplusHeNetHeating[i] =
- // (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
-
- // for (j = 0; j < cooling_N_He; j++) {
- // cooling_HplusHeNetHeating[i][j] =
- // (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
-
- // for (k = 0; k < cooling_N_nH; k++)
- // cooling_HplusHeNetHeating[i][j][k] =
- // (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
- // }
- //}
-
- float cooling_HplusHeNetHeating[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
-
- //cooling_HplusHeElectronAbundance =
- // (float ****)mymalloc("coolHHeNe", 2 * sizeof(float ***));
-
- //for (i = 0; i < 2; i++) {
- // cooling_HplusHeElectronAbundance[i] =
- // (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
-
- // for (j = 0; j < cooling_N_He; j++) {
- // cooling_HplusHeElectronAbundance[i][j] =
- // (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
-
- // for (k = 0; k < cooling_N_nH; k++)
- // cooling_HplusHeElectronAbundance[i][j][k] =
- // (float *)mymalloc("coolNH", cooling_N_Temp * sizeof(float));
- // }
- //}
-
- float cooling_HplusHeElectronAbundance[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
-
- //MakeNamePointers();
-
- /* Set H and He reionisation heating to erg/g */
- All.REION_H_Heating_ERGpG = All.REION_H_Heating_EVpH * EV_TO_ERG / PROTONMASS;
- All.REION_He_Heating_ERGpG =
- All.REION_He_Heating_EVpH * EV_TO_ERG / PROTONMASS;
-
- printf("Hydrogen reionisation heating = %e erg/g\n",
- All.REION_H_Heating_ERGpG);
- printf(" Helium reionisation heating = %e erg/g\n",
- All.REION_He_Heating_ERGpG);
-
- //cooling_ElementAbundance_SOLAR =
- // (double *)mymalloc("coolNele", cooling_N_Elements * sizeof(double));
- //cooling_ElementAbundance_SOLARM1 =
- // (double *)mymalloc("coolNele", cooling_N_Elements * sizeof(double));
-
- double cooling_ElementAbundance_SOLAR[cooling_N_Elements];
- double cooling_ElementAbundance_SOLARM1[cooling_N_Elements];
-
-#if EAGLE_TEST_COOL == 1
- if (ThisTask == 0)
- printf(
- "EAGLE_TEST_COOL option 1: ThisTask = 0 will evaluate cooling rates, "
- "then we will stop\n");
-
- DriverEagleTestCool();
-
- MPI_Barrier(MPI_COMM_WORLD);
- printf("Test 1 of cooling finished\n");
- endrun(-2);
-#endif
-
-#if EAGLE_TEST_COOL == 2
- if (ThisTask == 0)
- printf(
- "EAGLE_TEST_COOL option 2: ThisTask = 0 will evaluate thermal "
- "evolution at the mean density\n");
-
- DriverEagleTestCoolThermalEvolution();
-
- MPI_Barrier(MPI_COMM_WORLD);
- printf("Test 2 of cooling finished\n");
- endrun(-2);
-#endif
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine allocates header arrays for the cooling tables
- * ----------------------------------------------------------------------
- */
-
-void allocate_header_arrays(void) {
- int i;
-
- /* bins for interpolation in log10(T [K]) */
- cooling_Temp = (float *)mymalloc("coolTemp", cooling_N_Temp * sizeof(float));
-
- /* bins for interpolation in hydrogen number density, log10(nh [/cm^3]) */
- cooling_nH = (float *)mymalloc("coolNnH", cooling_N_nH * sizeof(float));
-
- /* bins for interpolation in thermal energy per unit mass, log10(u [erg/g]) */
- cooling_Therm = (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
-
- /* bins for interpolation in Helium abundance */
- cooling_HeFrac = (float *)mymalloc("coolNHe", cooling_N_He * sizeof(float));
-
- /* table of solar abundances */
- cooling_SolarAbundances =
- (float *)mymalloc("coolNSol", cooling_N_SolarAbundances * sizeof(float));
-
- /* assumed solar abundances used in constructing the tables, and corresponding
- * names */
- cooling_SolarAbundanceNames = (char **)mymalloc(
- "coolNsolName", cooling_N_SolarAbundances * sizeof(char *));
- if (ThisTask != 0)
- for (i = 0; i < cooling_N_SolarAbundances; i++)
- cooling_SolarAbundanceNames[i] =
- (char *)mymalloc("coolNsolName2", EL_NAME_LENGTH * sizeof(char));
-
- /* names of chemical elements stored in table */
- cooling_ElementNames =
- (char **)mymalloc("coolElNa", cooling_N_Elements * sizeof(char *));
- if (ThisTask != 0)
- for (i = 0; i < cooling_N_Elements; i++)
- cooling_ElementNames[i] =
- (char *)mymalloc("coolElNa2", EL_NAME_LENGTH * sizeof(char));
-}
-
-/*
- * ======================================================================
- * FUNCTIONS FOR BASIC TESTING OF COOLING ROUTINES - TEST_COOL flag.
- * ======================================================================
- */
-
-/*
- * ----------------------------------------------------------------------
- * This routine tests the interpolation of the cooling tables
- * ----------------------------------------------------------------------
- */
-
-#if EAGLE_TEST_COOL > 0
-void DriverEagleTestCool() {
- int iz;
-
- int z_index;
-
- float z1, z2, dz;
-
- /* we use redshifts at exact values from the tables, */
- /* as well as values in between */
- for (iz = cooling_N_Redshifts - 2; iz >= 0; iz--) {
- z1 = cooling_Redshifts[iz + 1]; /* collisional table */
- z2 = cooling_Redshifts[iz]; /* collisional table */
-
- get_redshift_index(z1, &z_index, &dz);
- LoadCoolingTables(z_index);
-
- /*
- EagleTestCool(z1 - 0.1*(z1-z2));
- */
-
- if (ThisTask == 0) EagleTestCool(z2);
- }
-}
-
-void EagleTestCool(float z) {
- FILE *outfile;
-
- int i, j, k, index1, index2, H_index, He_index, z_index;
-
- char filename[500];
-
- double Hefrac, ener, dens, temp;
-
- static double ParticleMetalFraction[EAGLE_NELEMENTS];
-
- float dz;
-
- /* given Helium abundance */
- Hefrac = cooling_HeFrac[0];
-
- get_redshift_index(z, &z_index, &dz);
-
- H_index = element_index("Hydrogen");
- He_index = element_index("Helium");
-
- /* loop over elements */
- for (i = 0; i < EAGLE_NELEMENTS; i++) {
- /* put abundance of element i to solar value, zero all others */
- for (j = 0; j < EAGLE_NELEMENTS; j++) ParticleMetalFraction[j] = 0;
-
- index1 = get_element_index(cooling_ElementNames, cooling_N_Elements,
- ElementNames[i]);
-
- if (index1 >= 0) /* index < 0 for Helium which is not a metal */
- {
- index2 = get_element_index(cooling_SolarAbundanceNames,
- cooling_N_SolarAbundances,
- cooling_ElementNames[index1]);
-
- if (index2 < 0 || index2 >= cooling_N_SolarAbundances) {
- printf(" sorry: element %s not found \n", cooling_ElementNames[index1]);
- endrun(-12344);
- }
-
- sprintf(filename, "testcool.z=%1.3f.Helium=%1.3f.Element=%s", z, Hefrac,
- ElementNames[i]);
- outfile = fopen(filename, "w");
- fprintf(outfile, "# element = %s\n", ElementNames[i]);
-
- ParticleMetalFraction[H_index] =
- 1.0 - Hefrac - cooling_SolarAbundances[index2];
- ParticleMetalFraction[He_index] = Hefrac;
-
- ParticleMetalFraction[i] = cooling_SolarAbundances[index2];
- } else {
- /* primordial cooling */
- sprintf(filename, "testcool.z=%1.3f.Helium=%1.3f.Element=%s", z, Hefrac,
- "Primordial");
- outfile = fopen(filename, "w");
- fprintf(outfile, "# element = %s\n", "Primordial");
-
- ParticleMetalFraction[H_index] = 1.0 - Hefrac;
- ParticleMetalFraction[He_index] = Hefrac;
- }
-
- /* write temperature density grid */
-
- fprintf(outfile, "# redshift = %e\n", z);
- fprintf(outfile, "# Helium abundance = %e\n", Hefrac);
- fprintf(outfile, "# cooling_N_Temp = %d\n", cooling_N_Temp);
- fprintf(outfile, "# cooling_N_nH = %d\n", cooling_N_nH);
- fprintf(outfile, "# T[K] u[erg/g] N_h [cm^-3] Lambda [erg cm^3 s^-1]\n");
- for (j = 0; j < cooling_N_Temp - 1; j++) {
- ener = cooling_Therm[j];
-
- for (k = 0; k < cooling_N_nH - 1; k++) {
- dens = cooling_nH[k];
-
- temp = eagle_convert_u_to_temp(dz, pow(10.0, ener), pow(10.0, dens),
- Hefrac);
-
- fprintf(outfile, " %e %e %e %e\n", temp, ener, pow(10.0, dens),
- EagleCoolingRate(pow(10.0, ener), pow(10.0, dens), z, dz,
- ParticleMetalFraction));
- }
- }
-
- fclose(outfile);
- }
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine tests the thermal evolution of gas at the critical density
- * ----------------------------------------------------------------------
- */
-
-void DriverEagleTestCoolThermalEvolution() {
- FILE *outfile = 0;
-
- char filename[100];
-
- int i, HeIndex, HIndex;
-
- int z_index;
-
- float z, dz;
-
- float zold, znow;
-
- double zmax = 49.0, deltaz = 0.1, uinit = 5.0e10, delta = 0.0, rho, u;
-
- double RhoCrit, hubble_z, dtime, da;
-
- double element_metallicity[EAGLE_NELEMENTS];
-
- double Hefrac = 0.248, Hfrac = 0.752, temp;
-
- double uc = 0;
-
- /* test only makes sense if co-moving integration is on */
- if (All.ComovingIntegrationOn != 1) {
- if (ThisTask == 0)
- printf(
- " this cooling test only works when All.ComovingIntegrationOn=1 \n");
- endrun(-2);
- }
-
- for (i = 0; i < EAGLE_NELEMENTS; i++) element_metallicity[i] = 0;
-
- HeIndex = element_index("Helium");
- HIndex = element_index("Hydrogen");
-
- element_metallicity[HeIndex] = Hefrac;
- element_metallicity[HIndex] = Hfrac;
-
- /* crit. density in code units */
- RhoCrit = 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G);
- /* crit. density in physical g/cm^3 */
- RhoCrit *= All.UnitDensity_in_cgs * All.HubbleParam * All.HubbleParam;
-
- u = uinit;
-
- if (ThisTask == 0) {
- sprintf(filename, "thermevol.data");
-
- if (!(outfile = fopen(filename, "w"))) {
- printf("can't open file %s\n", filename);
- endrun(888);
- }
-
- fprintf(outfile, "# delta = %g\n", delta);
- fprintf(outfile, "# uinit = %g\n", uinit);
-
- fprintf(outfile, "# metallicity = ");
-
- for (i = 0; i < EAGLE_NELEMENTS; i++)
- fprintf(outfile, " %g\t", element_metallicity[i]);
-
- fprintf(outfile, "\n");
- fprintf(outfile, "# z T[K] u[erg/g] rho[g/cm^3] reion heat [erg/g]\n");
- }
-
- for (z = zmax; z - deltaz > 0; z -= deltaz) {
- znow = z - deltaz;
- zold = z;
-
- if (ThisTask == 0) {
- printf(" REDSHIFT = %f, %f\n", zold, znow);
- fflush(stdout);
- }
-
- u *= pow(((1 + znow) / (1 + zold)), 2);
-
- get_redshift_index(znow, &z_index, &dz);
- LoadCoolingTables(z_index);
-
- /* physical baryon density at this z, for this over density, in g/cm^3 */
- rho = (1 + delta) * RhoCrit * All.OmegaBaryon * pow(1 + znow, 3);
-
- da = deltaz / pow(1 + znow, 2);
-
- hubble_z = (All.HubbleParam * HUBBLE) *
- sqrt(All.Omega0 * pow(1 + znow, 3) +
- (1 - All.OmegaLambda - All.Omega0) * pow(1 + znow, 2) +
- All.OmegaLambda);
-
- dtime = (1 + znow) * da / hubble_z;
-
- u = EagleDoCooling(u, rho, dtime, -deltaz, znow, dz, element_metallicity);
-
- /* extra energy from Hydrogen reionization */
- if (znow <= All.REION_H_ZCenter && znow + deltaz > All.REION_H_ZCenter) {
- u += All.REION_H_Heating_ERGpG;
- uc += All.REION_H_Heating_ERGpG;
- if (ThisTask == 0)
- printf(" Hydrogen reionzation: added %e erg/g at redshift %g\n", uc,
- znow);
- }
-
- /* extra energy from Helium reionization */
- /* Note that this term is added in the cooling calculation */
- /* u += eagle_helium_reionization_extraheat(znow, -deltaz); */
- uc += eagle_helium_reionization_extraheat(znow, -deltaz);
-
- temp = eagle_convert_u_to_temp(dz, u, rho * Hfrac / PROTONMASS, Hefrac);
-
- if (ThisTask == 0)
- fprintf(outfile, " %g \t %g \t %g \t %g \t %g\n", z, temp, u, rho, uc);
- }
-
- if (ThisTask == 0) fclose(outfile);
-}
-
-#endif /* EAGLE_TEST_COOL > 0 */
-
-#endif
diff --git a/src/cooling/EAGLE/eagle_interpol.h b/src/cooling/EAGLE/eagle_interpol.h
deleted file mode 100644
index 414c2467e030e7ca9015a59fcb5faa64fcac1d81..0000000000000000000000000000000000000000
--- a/src/cooling/EAGLE/eagle_interpol.h
+++ /dev/null
@@ -1,179 +0,0 @@
-#ifndef SWIFT_EAGLE_INTERPOL_H
-#define SWIFT_EAGLE_INTERPOL_H
-
-#include "cooling_read_table.h"
-
-/*
- * ----------------------------------------------------------------------
- * This routine returns the position i of a value x in a 1D table and the
- * displacement dx needed for the interpolation. The table is assumed to
- * be evenly spaced.
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE void get_index_1d(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);
- *dx = ((float)x - table[*i]) * dxm1;
- }
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE float interpol_1d(float *table, int i, float dx) {
- float result;
-
- result = (1 - dx) * table[i] + dx * table[i + 1];
-
- return result;
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table, int i, float dx) {
- double result;
-
- result = (1 - dx) * table[i] + dx * table[i + 1];
-
- return result;
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a bi-linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int j, float dx, float dy, int size) {
- float result;
- int index[4];
-
- index[0] = row_major_index_2d(i,j,size);
- index[1] = row_major_index_2d(i,j+1,size);
- index[2] = row_major_index_2d(i+1,j,size);
- index[3] = row_major_index_2d(i+1,j+1,size);
- 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]];
-
- return result;
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a bi-linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE double interpol_2d_dbl(double *table, int i, int j, double dx, double dy, int size) {
- double result;
- int index[4];
-
- index[0] = row_major_index_2d(i,j,size);
- index[1] = row_major_index_2d(i,j+1,size);
- index[2] = row_major_index_2d(i+1,j,size);
- index[3] = row_major_index_2d(i+1,j+1,size);
- 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]];
-
- return result;
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a tri-linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int j, int k, float dx, float dy,
- float dz, int size) {
- float result;
- int index[8];
-
- index[0] = row_major_index_3d(i,j,k,size);
- index[1] = row_major_index_3d(i,j,k+1,size);
- index[2] = row_major_index_3d(i,j+1,k,size);
- index[3] = row_major_index_3d(i,j+1,k+1,size);
- index[4] = row_major_index_3d(i+1,j,k,size);
- index[5] = row_major_index_3d(i+1,j,k+1,size);
- index[6] = row_major_index_3d(i+1,j+1,k,size);
- index[7] = row_major_index_3d(i+1,j+1,k+1,size);
- 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]] +
- (1 - dx) * dy * dz * table[index[3]] +
- dx * (1 - dy) * (1 - dz) * table[index[4]] +
- dx * (1 - dy) * dz * table[index[5]] +
- dx * dy * (1 - dz) * table[index[6]] +
- dx * dy * dz * table[index[7]];
-
- return result;
-}
-
-/*
- * ----------------------------------------------------------------------
- * This routine performs a quadri-linear interpolation
- * ----------------------------------------------------------------------
- */
-
-__attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int j, int k, int l, float dx,
- float dy, float dz, float dw, int size) {
- float result;
- int index[16];
-
- index[0] = row_major_index_4d(i,j,k,l,size);
- index[1] = row_major_index_4d(i,j,k,l+1,size);
- index[2] = row_major_index_4d(i,j,k+1,l,size);
- index[3] = row_major_index_4d(i,j,k+1,l+1,size);
- index[4] = row_major_index_4d(i,j+1,k,l,size);
- index[5] = row_major_index_4d(i,j+1,k,l+1,size);
- index[6] = row_major_index_4d(i,j+1,k+1,l,size);
- index[7] = row_major_index_4d(i,j+1,k+1,l+1,size);
- index[8] = row_major_index_4d(i+1,j,k,l,size);
- index[9] = row_major_index_4d(i+1,j,k,l+1,size);
- index[10] = row_major_index_4d(i+1,j,k+1,l,size);
- index[11] = row_major_index_4d(i+1,j,k+1,l+1,size);
- index[12] = row_major_index_4d(i+1,j+1,k,l,size);
- index[13] = row_major_index_4d(i+1,j+1,k,l+1,size);
- index[14] = row_major_index_4d(i+1,j+1,k+1,l,size);
- index[15] = row_major_index_4d(i+1,j+1,k+1,l+1,size);
-
- result = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table[index[0]] +
- (1 - dx) * (1 - dy) * (1 - dz) * dw * table[index[1]] +
- (1 - dx) * (1 - dy) * dz * (1 - dw) * table[index[2]] +
- (1 - dx) * (1 - dy) * dz * dw * table[index[3]] +
- (1 - dx) * dy * (1 - dz) * (1 - dw) * table[index[4]] +
- (1 - dx) * dy * (1 - dz) * dw * table[index[5]] +
- (1 - dx) * dy * dz * (1 - dw) * table[index[6]] +
- (1 - dx) * dy * dz * dw * table[index[7]] +
- dx * (1 - dy) * (1 - dz) * (1 - dw) * table[index[8]] +
- dx * (1 - dy) * (1 - dz) * dw * table[index[9]] +
- dx * (1 - dy) * dz * (1 - dw) * table[index[10]] +
- dx * (1 - dy) * dz * dw * table[index[11]] +
- dx * dy * (1 - dz) * (1 - dw) * table[index[12]] +
- dx * dy * (1 - dz) * dw * table[index[13]] +
- dx * dy * dz * (1 - dw) * table[index[14]] +
- dx * dy * dz * dw * table[index[15]];
-
- return result;
-}
-
-#endif
diff --git a/tests/testCooling.c b/tests/testCooling.c
index 7d077ef2e06d9208a13eb7ded71995a362606c92..be35f7969b907358665a3d3ca23691bed9e617c5 100644
--- a/tests/testCooling.c
+++ b/tests/testCooling.c
@@ -68,15 +68,24 @@ int main() {
exit(1);
}
+ gamma = 5.0/3.0;
+
+ chemistry_init(params, &us, &internal_const, &chemistry_data);
+ chemistry_first_init_part(&p,&xp,&chemistry_data);
+ chemistry_print(&chemistry_data);
+
+ u = 1.0*pow(10.0,11)/(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));
+
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
- chemistry_init(params, &us, &internal_const, &chemistry_data);
- chemistry_first_init_part(&p,&xp,&chemistry_data);
- chemistry_print(&chemistry_data);
for(int t_i = 0; t_i < n_t_i; t_i++){
@@ -86,8 +95,6 @@ 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));
- 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",temperature_cgs,cooling_du_dt);