diff --git a/examples/CoolingBox/README b/examples/CoolingBox/README
new file mode 100644
index 0000000000000000000000000000000000000000..1e41bad4afaf690f41a8c75f05b0452fea51b022
--- /dev/null
+++ b/examples/CoolingBox/README
@@ -0,0 +1,5 @@
+This is a test of cooling on a box of uniform gas. It can be run with a variety of of cooling models, both with and without cosmology. If running without cosmology specify coolingBox_non-cosmo.yml
+
+In the case of EAGLE cooling a script (test_script.sh) is provided to run with a variety of parameters and the solution verified against a subcycled explicit solution calculated by the python script analytical_test.py. The parameters that are tested for are: redshift, metal abundance (solar or primordial), hydrogen number density and initial energy. These parameters may be changed in test_script.sh by modifying the for loops, or specifying the desired metal abundances. Currently the script is designed to work with cosmology only. Please note that in order for the cooling box solution to be verified against the subcycled solution, the cooling rates for the parameters specified in test_script.sh are calculated by ../CoolingRates/testCooling.
+
+Initial conditions are produced from a glass file, specified in run.sh and makeIC.py.
diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 16c2fc55777eafb367896b4454f2c716538e6058..23387bcf91b294230ef4ce8b01d117e706e0f7af 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/examples/CoolingBox/analytical_test.py
@@ -136,15 +136,15 @@ log_nh = float(sys.argv[2])
redshift = float(sys.argv[1])
p = plt.figure(figsize = (6,6))
p1, = plt.loglog(t_sol,u_sol,linewidth = 0.5,color = 'k', marker = '*', markersize = 1.5,label = 'explicit ODE')
-p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean alexei')
+p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean')
l = legend(handles = [p1,p2])
xlabel("${\\rm{Time~[s]}}$", labelpad=0, fontsize = 14)
ylabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
if int(sys.argv[4]) == 1:
- title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error alexei: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
+ title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_solar.png"
elif int(sys.argv[4]) == 0:
- title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error alexei: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
+ title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_zero_metal.png"
savefig(name, dpi=200)
diff --git a/examples/CoolingBox/coolingBox_non-cosmo.yml b/examples/CoolingBox/coolingBox_non-cosmo.yml
new file mode 100644
index 0000000000000000000000000000000000000000..d8627e8bd2072bfcb20e7abd5aedda7b336566cb
--- /dev/null
+++ b/examples/CoolingBox/coolingBox_non-cosmo.yml
@@ -0,0 +1,80 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 2.0e33 # Solar masses
+ UnitLength_in_cgs: 3.0857e21 # Kiloparsecs
+ UnitVelocity_in_cgs: 1.0e5 # Kilometers per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.01 # The end time of the simulation (in internal units).
+ dt_min: 1e-9 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-5 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: coolingBox # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 2e-5 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./coolingBox.hdf5 # The file to read
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+ lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
+ minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
+ mean_molecular_weight: 0.59 # Mean molecular weight
+ hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
+ cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
+
+# Cooling with Grackle 2.0
+GrackleCooling:
+ CloudyTable: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ WithUVbackground: 0 # Enable or not the UV background
+ Redshift: 0 # Redshift to use (-1 means time based redshift)
+ WithMetalCooling: 1 # Enable or not the metal cooling
+ ProvideVolumetricHeatingRates: 0 # User provide volumetric heating rates
+ ProvideSpecificHeatingRates: 0 # User provide specific heating rates
+ SelfShieldingMethod: 0 # Grackle (<= 3) or Gear self shielding method
+ OutputMode: 1 # Write in output corresponding primordial chemistry mode
+ MaxSteps: 1000
+ ConvergenceLimit: 1e-2
+
+EagleCooling:
+ filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
+ reionisation_redshift: 8.898
+ he_reion: 1
+ he_reion_z_center: 3.5
+ he_reion_z_sigma: 0.5
+ he_reion_ev_pH: 2.0
+
+EAGLEChemistry:
+ InitMetallicity: 0.014
+ InitAbundance_Hydrogen: 0.70649785
+ InitAbundance_Helium: 0.28055534
+ InitAbundance_Carbon: 2.0665436e-3
+ InitAbundance_Nitrogen: 8.3562563e-4
+ InitAbundance_Oxygen: 5.4926244e-3
+ InitAbundance_Neon: 1.4144605e-3
+ InitAbundance_Magnesium: 5.907064e-4
+ InitAbundance_Silicon: 6.825874e-4
+ InitAbundance_Iron: 1.1032152e-3
+ CalciumOverSilicon: 0.0941736
+ SulphurOverSilicon: 0.6054160
+
+GearChemistry:
+ InitialMetallicity: 0.01295
+
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 03f5626023bf1570caa65df52673c8b5565ec4a3..eb76f3d7b5b7a100fdcaf9b38b5d4c82f2f20104 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -27,8 +27,8 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec
-rho = 34504.4797588 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 44554455.4455 # Pressure in code units (at 10^6K)
+rho = 3450447.97588 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+P = 4455445544.55 # 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"
diff --git a/examples/CoolingRates/README b/examples/CoolingRates/README
new file mode 100644
index 0000000000000000000000000000000000000000..414a4dc9ebb453cd0a1c59160a36139b437c025f
--- /dev/null
+++ b/examples/CoolingRates/README
@@ -0,0 +1,7 @@
+This is a test that produces a plot of the contribution to the cooling rate from each of the elements depending on internal energy, density and redshift based on the EAGLE tables. To do so, the function in src/cooling/EAGLE returning the cooling rate is run for multiple values of the internal energy. The resulting cooling rates are written to files and plotted with a python script (cooling_rates_plot.py).
+
+The test may be run by using the shell script, test_script.sh, or by running:
+./testCooling -z X -d Y
+python cooling_rates_plot.py
+
+where X is the redshift at which the cooling rates are evaluated and Y is the base 10 logarithm of the hydrogen number density. Different metallicities may be specified in testCooling.yml
diff --git a/examples/CoolingRates/cooling_rates_plot.py b/examples/CoolingRates/cooling_rates_plot.py
new file mode 100644
index 0000000000000000000000000000000000000000..d115ad5e580439f1e8e65338a19f717497f34d34
--- /dev/null
+++ b/examples/CoolingRates/cooling_rates_plot.py
@@ -0,0 +1,39 @@
+import matplotlib.pyplot as plt
+
+elements = 11
+u = []
+cooling_rate = [[] for i in range(elements+1)]
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+ data = line.split()
+ u.append(float(data[0]))
+ cooling_rate[0].append(-float(data[1]))
+
+file_in.close()
+
+for elem in range(elements):
+ file_in = open('cooling_element_'+str(elem)+'.dat','r')
+ for line in file_in:
+ data = line.split()
+ cooling_rate[elem+1].append(-float(data[0]))
+ file_in.close()
+
+ax = plt.subplot(111)
+p0, = plt.loglog(u, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(u, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
+p2, = plt.loglog(u, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
+p3, = plt.loglog(u, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p4, = plt.loglog(u, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p5, = plt.loglog(u, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
+p6, = plt.loglog(u, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p7, = plt.loglog(u, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
+p8, = plt.loglog(u, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p9, = plt.loglog(u, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p10, = plt.loglog(u, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
+ax.set_position([0.15,0.15,0.75,0.75])
+plt.xlim([1e12,1e17])
+plt.ylim([1e-24,1e-21])
+plt.xlabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
+plt.ylabel("${\Lambda/n_h^2 }$ ${\\rm{[erg \cdot cm^3 \cdot s^{-1}]}}$", fontsize = 14)
+plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
+plt.show()
diff --git a/examples/CoolingRates/testCooling.c b/examples/CoolingRates/testCooling.c
index 4f141097a99d1b73e4eac4c6ffb3f3098cebdd9f..22ca23423ebc87cf8ca5a9ba0b08762641df13ab 100644
--- a/examples/CoolingRates/testCooling.c
+++ b/examples/CoolingRates/testCooling.c
@@ -69,299 +69,6 @@ void set_quantities(struct part *restrict p,
if(cosmo->z >= 1) hydro_set_init_internal_energy(p,(u*pow(scale_factor,2))/cooling->internal_energy_scale);
}
-/*
- * @brief Construct 1d tables from 4d EAGLE tables by
- * interpolating over redshift, hydrogen number density
- * and helium fraction.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param temp_table Pointer to 1d interpolated table of temperature values
- * @param H_plus_He_heat_table Pointer to 1d interpolated table of cooling rates
- * due to hydrogen and helium
- * @param H_plus_He_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to hydrogen and helium
- * @param element_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to metals
- * @param element_print_cooling_table Pointer to 1d interpolated table of
- * cooling rates due to each of the metals
- * @param element_cooling_table Pointer to 1d interpolated table of cooling
- * rates due to the contribution of all the metals
- * @param abundance_ratio Pointer to array of ratios of metal abundances to solar
- * @param ub Upper bound in temperature on table construction
- * @param lb Lower bound in temperature on table construction
- */
-void construct_1d_tables_test(const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict internal_const,
- double *temp_table,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_electron_abundance_table,
- double *element_print_cooling_table,
- double *element_cooling_table,
- float *abundance_ratio,
- float *ub, float *lb){
-
- // obtain mass fractions and number density for particle
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, helium fraction, hydrogen number density indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- if (cosmo->z > cooling->reionisation_redshift) {
- // Photodissociation table
- printf("Eagle testCooling.c photodissociation table redshift reionisation redshift %.5e %.5e\n", cosmo->z, cooling->reionisation_redshift);
- construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.temperature,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,
- He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.H_plus_He_electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
- cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_2d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
- printf("Eagle testCooling.c no compton table redshift max redshift %.5e %.5e\n", cosmo->z, cooling->Redshifts[cooling->N_Redshifts - 1]);
- // High redshift table
- construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.temperature,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,
- He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.H_plus_He_electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
- cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_2d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- } else {
- // Normal tables
- printf("Eagle testCooling.c normal table redshift %.5e\n", cosmo->z);
- construct_1d_table_from_4d(p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.temperature,
- z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
- cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- }
-
-}
-
-/*
- * @brief Compare calculation of dlambda/du (gradient of cooling rate,
- * needed for Newton's method) between interpolating 1d tables and
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_dlambda_du(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- // allocate tables
- double H_plus_He_heat_table[cooling->N_Temp];
- double H_plus_He_electron_abundance_table[cooling->N_Temp];
- double temp_table[cooling->N_Temp];
- double element_cooling_table[cooling->N_Temp];
- double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];
- double element_electron_abundance_table[cooling->N_Temp];
- double rate_element_table[cooling->N_Elements+2];
- float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
- double dlambda_du1, dlambda_du2;
-
- // calculate ratio of particle metal abundances to solar
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- // set hydrogen number density and internal energy
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
- // extract hydrogen and helium mass fractions
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, hydrogen number density and helium fraction indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-
- // construct tables
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling->Temp[0]/eagle_log_10_e;
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-
- // calculate dlambda/du for different values of internal energy
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- p,cooling,cosmo,internal_const,rate_element_table);
- cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
- &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- p,cooling,cosmo,internal_const,NULL,
- abundance_ratio);
- printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
- }
-}
-
-/*
- * @brief Compare calculating temperature between interpolating 1d tables and
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_temp(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- // allocate tables
- double H_plus_He_heat_table[cooling->N_Temp];
- double H_plus_He_electron_abundance_table[cooling->N_Temp];
- double temp_table[cooling->N_Temp];
- double element_cooling_table[cooling->N_Temp];
- double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];
- double element_electron_abundance_table[cooling->N_Temp];
- float *abundance_ratio;
-
- // calculate ratio of particle metal abundances to solar
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- // set hydrogen number density and internal energy
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
- // extract hydrogen and helium mass fractions
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, hydrogen number density and helium fraction indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- // construct tables
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling->Temp[0]/eagle_log_10_e;
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-
- // calculate temperature for different values of internal energy
- float T1d, T4d, delta_u;
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- cooling);
- T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He,cooling, cosmo);
- printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
- }
-}
-
/*
* @brief Produces contributions to cooling rates for different
* hydrogen number densities, from different metals,
@@ -447,14 +154,6 @@ int main(int argc, char **argv) {
abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
- // Declare 1D tables
- double H_plus_He_heat_table[cooling.N_Temp];
- double H_plus_He_electron_abundance_table[cooling.N_Temp];
- double temp_table[cooling.N_Temp];
- double element_cooling_table[cooling.N_Temp];
- double element_print_cooling_table[cooling.N_Elements * cooling.N_Temp];
- double element_electron_abundance_table[cooling.N_Temp];
-
// extract mass fractions, calculate table indices and offsets
float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
@@ -464,139 +163,48 @@ int main(int argc, char **argv) {
get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
- float upper_bound = cooling.Temp[cooling.N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling.Temp[0]/eagle_log_10_e;
-
- int nt = 250;//, n_nh = 6;
+ int nt = 250;
double u = pow(10.0,10);
- //if (argc == 1 || strcmp(argv[1], "nh") == 0){
- // Calculate cooling rates at different densities
- //for(int i = 0; i < n_nh; i++){
- // // Open files
- // char output_filename[21];
- // sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
- // FILE *output_file = fopen(output_filename, "w");
- // if (output_file == NULL) {
- // printf("Error opening file!\n");
- // exit(1);
- // }
-
- // // set hydrogen number density, construct tables
- // nh = pow(10.0,-i);
- // set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- // construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- // temp_table, H_plus_He_heat_table,
- // H_plus_He_electron_abundance_table,
- // element_electron_abundance_table,
- // element_print_cooling_table,
- // element_cooling_table,
- // abundance_ratio, &upper_bound, &lower_bound);
- // get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
- // for(int j = 0; j < nt; j++){
- // // set internal energy
- // set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
- // u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- // double dlambda_du;
- // float delta_u, cooling_du_dt, logT;
-
- // // calculate cooling rates using 1d tables
- // if (tables == 1) {
- // cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
- // H_plus_He_electron_abundance_table,
- // element_print_cooling_table,
- // element_electron_abundance_table,
- // temp_table,
- // z_index, dz, n_h_i, d_n_h, He_i, d_He,
- // &p,&cooling,&cosmo,&internal_const,NULL,
- // abundance_ratio);
- // logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- // &p,&cooling, &cosmo, &internal_const);
- // } else {
- // // calculate cooling rates using 4d tables
- // cooling_du_dt = eagle_metal_cooling_rate(log10(u),
- // &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- // &p,&cooling,&cosmo,&internal_const,NULL,
- // abundance_ratio);
- // logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- // dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
- // }
- // float temperature_swift = pow(10.0,logT);
- // fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- // }
- // fclose(output_file);
- //}
- //}
// Calculate contributions from metals to cooling rate
- //if (argc >= 1 || strcmp(argv[1],"metals") == 0) {
- // open file
- char output_filename[21];
- sprintf(output_filename, "%s", "cooling_output.dat");
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density, construct 1d tables
- if (log_10_nh == 100) {
- nh = 1.0e-1;
- } else {
- nh = pow(10.0,log_10_nh);
- }
- //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
- printf("Eagle testcooling.c nh %.5e\n", nh);
- u = pow(10.0,14.0);
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- temp_table, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_electron_abundance_table,
- element_print_cooling_table,
- element_cooling_table,
- abundance_ratio, &upper_bound, &lower_bound);
- float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- printf("Eagle testcooling.c inn_h %.5e\n", inn_h);
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- // Loop over internal energy
- for(int j = 0; j < nt; j++){
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
- u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- float delta_u, cooling_du_dt, logT;
-
- // calculate cooling rates using 1d tables
- if (tables == 1) {
- cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- &p,&cooling,&cosmo,&internal_const);
- logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- &cooling);
- } else {
- // calculate cooling rates using 4d tables
- cooling_du_dt = eagle_print_metal_cooling_rate(
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,
- abundance_ratio);
- logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He,&cooling, &cosmo);
- }
- //float temperature_swift = pow(10.0,logT);
- //fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
- }
- fclose(output_file);
- //}
+ // open file
+ char output_filename[21];
+ sprintf(output_filename, "%s", "cooling_output.dat");
+ FILE *output_file = fopen(output_filename, "w");
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
- // compare temperatures and dlambda/du calculated from 1d and 4d tables
- //compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
- //compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
+ // set hydrogen number density, construct 1d tables
+ if (log_10_nh == 100) {
+ nh = 1.0e-1;
+ } else {
+ nh = pow(10.0,log_10_nh);
+ }
+ //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
+ u = pow(10.0,14.0);
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
+ float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+ &internal_const) *
+ cooling.number_density_scale;
+ get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
+
+ // Loop over internal energy
+ for(int j = 0; j < nt; j++){
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
+ u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
+ float cooling_du_dt;
+
+ // calculate cooling rates using 4d tables
+ cooling_du_dt = eagle_print_metal_cooling_rate(
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,
+ abundance_ratio);
+ fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
+ }
+ fclose(output_file);
+ printf("done\n");
free(params);
return 0;
diff --git a/examples/CoolingRates/testCooling.yml b/examples/CoolingRates/testCooling.yml
index db3a7cc00a529dcc403d4c24c11b30186baa7896..3010636433d35421009539d1858722ea89820bc9 100644
--- a/examples/CoolingRates/testCooling.yml
+++ b/examples/CoolingRates/testCooling.yml
@@ -59,16 +59,16 @@ InitialConditions:
cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
EAGLEChemistry:
- InitMetallicity: 0.0
- InitAbundance_Hydrogen: 0.752
- InitAbundance_Helium: 0.248
- InitAbundance_Carbon: 0.000
- InitAbundance_Nitrogen: 0.000
- InitAbundance_Oxygen: 0.000
- InitAbundance_Neon: 0.000
- InitAbundance_Magnesium: 0.000
- InitAbundance_Silicon: 0.000
- InitAbundance_Iron: 0.000
+ InitMetallicity: 0.014
+ InitAbundance_Hydrogen: 0.70649785
+ InitAbundance_Helium: 0.28055534
+ InitAbundance_Carbon: 2.0665436e-3
+ InitAbundance_Nitrogen: 8.3562563e-4
+ InitAbundance_Oxygen: 5.4926244e-3
+ InitAbundance_Neon: 1.4144605e-3
+ InitAbundance_Magnesium: 5.907064e-4
+ InitAbundance_Silicon: 6.825874e-4
+ InitAbundance_Iron: 1.1032152e-3
CalciumOverSilicon: 0.0941736
SulphurOverSilicon: 0.6054160
diff --git a/examples/CoolingRates/test_script.sh b/examples/CoolingRates/test_script.sh
new file mode 100644
index 0000000000000000000000000000000000000000..85847710da33db97d30b9e3aa7f361f44f18f55a
--- /dev/null
+++ b/examples/CoolingRates/test_script.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+./testCooling -z 0.1 -d -3
+python cooling_rates_plot.py
diff --git a/examples/CoolingTest_Alexei/Makefile.am b/examples/CoolingTest_Alexei/Makefile.am
deleted file mode 100644
index e878a3792a84f7ace5f3b60fcc5875169145731b..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/Makefile.am
+++ /dev/null
@@ -1,32 +0,0 @@
-# tHIS FIle is part of SWIFT.
-# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
-#
-# This program is free software: you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation, either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program. If not, see <http://www.gnu.org/licenses/>.
-
-# Add the source directory and the non-standard paths to the included library headers to CFLAGS
-AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
-
-AM_LDFLAGS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
-
-# Extra libraries.
-EXTRA_LIBS = $(HDF5_LIBS)
-
-# Programs.
-bin_PROGRAMS = testCooling
-
-# Sources
-testCooling_SOURCES = testCooling.c
-testCooling_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
-testCooling_LDADD = ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
-
diff --git a/examples/CoolingTest_Alexei/plots.py b/examples/CoolingTest_Alexei/plots.py
deleted file mode 100644
index 969e6de713a4c421cc14d5f88ef68bafaab59a41..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/plots.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 11
-temperature = []
-cooling_rate = [[] for i in range(elements+1)]
-u = []
-fu = []
-length = 0
-file_in = open('cooling_output.dat', 'r')
-for line in file_in:
- data = line.split()
- temperature.append(float(data[0]))
- cooling_rate[0].append(-float(data[1]))
-
-file_in.close()
-
-for elem in range(elements):
- file_in = open('cooling_element_'+str(elem)+'.dat','r')
- for line in file_in:
- data = line.split()
- cooling_rate[elem+1].append(-float(data[0]))
- file_in.close()
-
-#file_in = open('newton_output.dat', 'r')
-#for line in file_in:
-# data = line.split()
-# u.append(float(data[0]))
-# fu.append(float(data[1]))
-#
-#file_in.close()
-
-#p0, = plt.plot(u,fu, color = 'k')
-#p1, = plt.plot(u,[0 for x in u], color = 'r')
-#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
-#plt.xlim([1e13,2.0e14])
-#plt.ylim([0e13,2e14])
-#plt.xlabel('u')
-#plt.ylabel('f(u)')
-#plt.show()
-
-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 = '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([1e2,1e8])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
-#plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
-plt.show()
diff --git a/examples/CoolingTest_Alexei/plots_nh.py b/examples/CoolingTest_Alexei/plots_nh.py
deleted file mode 100644
index ab3300750e2497ec8026dbe2d5c6c7d358f87831..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/plots_nh.py
+++ /dev/null
@@ -1,47 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 6
-temperature = [[] for i in range(elements+1)]
-cooling_rate = [[] for i in range(elements+1)]
-u = []
-fu = []
-length = 0
-
-for elem in range(elements):
- file_in = open('cooling_output_'+str(elem)+'.dat','r')
- for line in file_in:
- data = line.split()
- temperature[elem+1].append(float(data[0]))
- cooling_rate[elem+1].append(-float(data[1]))
- file_in.close()
-
-#file_in = open('newton_output.dat', 'r')
-#for line in file_in:
-# data = line.split()
-# u.append(float(data[0]))
-# fu.append(float(data[1]))
-#
-#file_in.close()
-
-#p0, = plt.plot(u,fu, color = 'k')
-#p1, = plt.plot(u,[0 for x in u], color = 'r')
-#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
-#plt.xlim([1e13,2.0e14])
-#plt.ylim([0e13,2e14])
-#plt.xlabel('u')
-#plt.ylabel('f(u)')
-#plt.show()
-
-#p0, = plt.loglog(temperature[0], cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
-p1, = plt.loglog(temperature[1], cooling_rate[1], linewidth = 0.5, color = 'k', label = 'nh = 10^0')
-p2, = plt.loglog(temperature[2], cooling_rate[2], linewidth = 0.5, color = 'b', label = 'nh = 10^-1')
-p3, = plt.loglog(temperature[3], cooling_rate[3], linewidth = 0.5, color = 'g', label = 'nh = 10^-2')
-p4, = plt.loglog(temperature[4], cooling_rate[4], linewidth = 0.5, color = 'r', label = 'nh = 10^-3')
-p5, = plt.loglog(temperature[5], cooling_rate[5], linewidth = 0.5, color = 'c', label = 'nh = 10^-4')
-p6, = plt.loglog(temperature[6], cooling_rate[6], linewidth = 0.5, color = 'm', label = 'nh = 10^-5')
-plt.ylim([1e-24,1e-21])
-plt.xlim([1e4,1e8])
-plt.legend(handles = [p1,p2,p3,p4,p5,p6])
-plt.xlabel('Temperature (K)')
-plt.ylabel('Cooling rate (eagle units...)')
-plt.show()
diff --git a/examples/CoolingTest_Alexei/testCooling.c b/examples/CoolingTest_Alexei/testCooling.c
deleted file mode 100644
index 3a8cc301812097baf1d9e136887de86000097198..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling.c
+++ /dev/null
@@ -1,420 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-enum table_index {
- EAGLE_Hydrogen = 0,
- EAGLE_Helium,
- EAGLE_Carbon,
- EAGLE_Nitrogen,
- EAGLE_Oxygen,
- EAGLE_Neon,
- EAGLE_Magnesium,
- EAGLE_Silicon,
- EAGLE_Iron
-};
-
-void set_quantities(struct part *restrict p,
- const struct unit_system *restrict us,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict internal_const,
- float nh,
- double u){
-
- const float gamma = 5.0/3.0;
- double hydrogen_number_density = nh*pow(units_cgs_conversion_factor(us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo->z,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]);
-
- float scale_factor = 1.0/(1.0+cosmo->z);
- float pressure = (u*pow(scale_factor,3))/cooling->internal_energy_scale*p->rho*(gamma -1.0);
- p->entropy = pressure/(pow(p->rho,gamma));
-}
-
-void construct_1d_tables_test(const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict internal_const,
- double *temp_table,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_electron_abundance_table,
- double *element_print_cooling_table,
- double *element_cooling_table,
- float *abundance_ratio,
- float *ub, float *lb){
-
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
-
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
- printf("testCooling 1d z_index dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n", z_index,dz, n_h_i, d_n_h, He_i, d_He);
-
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.temperature,
- z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
- cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table,abundance_ratio, ub, lb);
- construct_1d_print_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table,abundance_ratio, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-
-}
-
-void compare_dlambda_du(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double temp_table[176];
- double element_cooling_table[176];
- double element_print_cooling_table[9 * 176];
- double element_electron_abundance_table[176];
- double rate_element_table[11];
- float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
- double dlambda_du1, dlambda_du2;
-
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
- const float log_10_e = 0.43429448190325182765;
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
- float lower_bound = cooling->Temp[0]/log_10_e;
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
- printf("testCooling.c 4d z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- p,cooling,cosmo,internal_const,rate_element_table,
- abundance_ratio);
- cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
- &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- p,cooling,cosmo,internal_const,NULL,
- abundance_ratio);
- printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
- }
-}
-
-void compare_temp(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double temp_table[176];
- double element_cooling_table[176];
- double element_print_cooling_table[9 * 176];
- double element_electron_abundance_table[176];
- float *abundance_ratio;
-
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
-
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
-
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- const float log_10_e = 0.43429448190325182765;
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
- float lower_bound = cooling->Temp[0]/log_10_e;
-
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
- float T1d, T4d, delta_u;
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
- printf("testCooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- p,cooling, cosmo, internal_const);
- T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, p,cooling, cosmo, internal_const);
- printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
- }
-}
-
-int main(int argc, char **argv) {
- /* Declare relevant structs */
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const internal_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./testCooling.yml";
-
- int tables = 0;
-
- /* Read the parameter file */
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- /* And dump the parameters as used. */
- parser_write_params_to_file(params, "used_parameters.yml");
-
- /* Init units */
- units_init_from_params(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &internal_const);
-
- /* Init chemistry */
- chemistry_init(params, &us, &internal_const, &chem_data);
- chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- /* Init cosmology */
- cosmology_init(params, &us, &internal_const, &cosmo);
- cosmology_print(&cosmo);
- cosmo.z = 0.9090909;
-
- /* Init cooling */
- cooling_init(params, &us, &internal_const, &cooling);
- cooling_print(&cooling);
-
- /* Calculate abundance ratios */
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- /* Declare 1D tables */
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double temp_table[176];
- double element_cooling_table[176];
- double element_print_cooling_table[9 * 176];
- double element_electron_abundance_table[176];
-
- float nh;
-
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
- const float log_10_e = 0.43429448190325182765;
- float upper_bound = cooling.Temp[cooling.N_Temp-1]/log_10_e;
- float lower_bound = cooling.Temp[0]/log_10_e;
-
- /* Loop over densities */
- int nt = 250, n_nh = 6;
- double u = pow(10.0,10);
- if (argc == 1 || strcmp(argv[1], "nh") == 0){
- for(int i = 0; i < n_nh; i++){
- /* Open files */
- char output_filename[21];
- sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- nh = pow(10.0,-i);
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- temp_table, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_electron_abundance_table,
- element_print_cooling_table,
- element_cooling_table,
- abundance_ratio, &upper_bound, &lower_bound);
- get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
- for(int j = 0; j < nt; j++){
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
- u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- double dlambda_du;
- float delta_u, cooling_du_dt, logT;
- if (tables == 1) {
- cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,NULL,
- abundance_ratio);
- logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- &p,&cooling, &cosmo, &internal_const);
- } else {
- cooling_du_dt = eagle_metal_cooling_rate(log10(u),
- &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,NULL,
- abundance_ratio);
- logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
- }
- float temperature_swift = pow(10.0,logT);
-
- fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- }
- fclose(output_file);
- }
- }
- if (argc == 1 || strcmp(argv[1],"metals") == 0) {
- printf("here \n");
- char output_filename[21];
- sprintf(output_filename, "%s", "cooling_output.dat");
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density, construct 1d tables
- //nh = pow(10.0,strtod(argv[2],NULL));
- nh = 1.0e-1;
- u = pow(10.0,14.0);
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- float nh_swift = chemistry_get_number_density(&p, &cosmo, chemistry_element_H, &internal_const)*cooling.number_density_scale;
- printf("hydrogen number density %.5e\n", nh_swift);
- construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- temp_table, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_electron_abundance_table,
- element_print_cooling_table,
- element_cooling_table,
- abundance_ratio, &upper_bound, &lower_bound);
- float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
- printf("testcooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n",z_index, dz, n_h_i, d_n_h, He_i, d_He);
- for(int j = 0; j < nt; j++){
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
- u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- float delta_u, cooling_du_dt, logT;
- if (tables == 1) {
- cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,
- abundance_ratio);
- logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- &p,&cooling, &cosmo, &internal_const);
- } else {
- cooling_du_dt = eagle_print_metal_cooling_rate(
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,
- abundance_ratio);
- }
- fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
- }
- fclose(output_file);
- }
-
- if (argc == 1 || strcmp(argv[1],"compare_temp") == 0) compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
- if (argc == 1 || strcmp(argv[1],"compare_dlambda_du") == 0) compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-
- free(params);
- return 0;
-}
-
-
diff --git a/examples/CoolingTest_Alexei/testCooling.yml b/examples/CoolingTest_Alexei/testCooling.yml
deleted file mode 100644
index f795ed6836786ae1ad72234fda8184c946e4ebe9..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling.yml
+++ /dev/null
@@ -1,94 +0,0 @@
-# Define the system of units to use internally.
-InternalUnitSystem:
- UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
- UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
- UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
- UnitCurrent_in_cgs: 1 # Amperes
- UnitTemp_in_cgs: 1 # Kelvin
-
-# Cosmological parameters
-Cosmology:
- h: 0.6777 # Reduced Hubble constant
- #a_begin: 0.9090909 # Initial scale-factor of the simulation
- a_begin: 1.0 # Initial scale-factor of the simulation
- a_end: 1.0 # Final scale factor of the simulation
- Omega_m: 0.307 # Matter density parameter
- Omega_lambda: 0.693 # Dark-energy density parameter
- Omega_b: 0.0455 # Baryon density parameter
-
-# Parameters governing the time integration
-TimeIntegration:
- time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1e-2 # The end time of the simulation (in internal units).
- dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
-
-Scheduler:
- max_top_level_cells: 15
-
-# Parameters governing the snapshots
-Snapshots:
- basename: eagle # Common part of the name of output files
- scale_factor_first: 0.92 # Scale-factor of the first snaphot (cosmological run)
- time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
- delta_time: 1.10 # Time difference between consecutive outputs (in internal units)
- compression: 1
-
-# Parameters governing the conserved quantities statistics
-Statistics:
- scale_factor_first: 0.92 # Scale-factor of the first stat dump (cosmological run)
- time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
- delta_time: 1.05 # Time between statistics output
-
-# Parameters for the self-gravity scheme
-Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
- comoving_softening: 0.0026994 # Comoving softening length (in internal units).
- max_physical_softening: 0.0007 # Physical softening length (in internal units).
-
-# Parameters for the hydrodynamics scheme
-SPH:
- resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
- CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
- minimal_temperature: 100 # (internal units)
-
-# Parameters related to the initial conditions
-InitialConditions:
- file_name: ./EAGLE_ICs_12.hdf5 # The file to read
- cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
-
-EAGLEChemistry:
- InitMetallicity: 0.014
- InitAbundance_Hydrogen: 0.70649785
- InitAbundance_Helium: 0.28055534
- InitAbundance_Carbon: 2.0665436e-3
- InitAbundance_Nitrogen: 8.3562563e-4
- InitAbundance_Oxygen: 5.4926244e-3
- InitAbundance_Neon: 1.4144605e-3
- InitAbundance_Magnesium: 5.907064e-4
- InitAbundance_Silicon: 6.825874e-4
- InitAbundance_Iron: 1.1032152e-3
- CalciumOverSilicon: 0.0941736
- SulphurOverSilicon: 0.6054160
- #InitMetallicity: 0.
- #InitAbundance_Hydrogen: 0.752
- #InitAbundance_Helium: 0.248
- #InitAbundance_Carbon: 0.000
- #InitAbundance_Nitrogen: 0.000
- #InitAbundance_Oxygen: 0.000
- #InitAbundance_Neon: 0.000
- #InitAbundance_Magnesium: 0.000
- #InitAbundance_Silicon: 0.000
- #InitAbundance_Iron: 0.000
- #CalciumOverSilicon: 0.0941736
- #SulphurOverSilicon: 0.6054160
-
-EagleCooling:
- filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
- reionisation_redshift: 20
- he_reion: 0
- he_reion_z_center: 3.5
- he_reion_z_sigma: 0.5
- he_reion_ev_pH: 2.0
-
diff --git a/examples/CoolingTest_Alexei/testCooling_messy.c b/examples/CoolingTest_Alexei/testCooling_messy.c
deleted file mode 100644
index b5c00e5625c43e57e46c5f8426640b9f511bc58a..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling_messy.c
+++ /dev/null
@@ -1,354 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main(int argc, char **argv) {
-
- /* Read the parameter file */
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const internal_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./testCooling.yml";
- enum table_index {
- EAGLE_Hydrogen = 0,
- EAGLE_Helium,
- EAGLE_Carbon,
- EAGLE_Nitrogen,
- EAGLE_Oxygen,
- EAGLE_Neon,
- EAGLE_Magnesium,
- EAGLE_Silicon,
- EAGLE_Iron
- };
-
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- /* And dump the parameters as used. */
- parser_write_params_to_file(params, "used_parameters.yml");
-
- double UnitMass_in_cgs = 1.989e43;
- double UnitLength_in_cgs = 3.085678e24;
- double UnitVelocity_in_cgs = 1e5;
- double UnitCurrent_in_cgs = 1;
- double UnitTemp_in_cgs = 1;
- units_init(&us, UnitMass_in_cgs, UnitLength_in_cgs, UnitVelocity_in_cgs, UnitCurrent_in_cgs, UnitTemp_in_cgs);
- phys_const_init(&us, params, &internal_const);
-
- double hydrogen_number_density_cgs = 1e-4;
- double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
- temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
- u_cgs = 0;
- cooling_du_dt = 0;
- temperature_cgs = 0;
- newton_func = 0;
- // int n_t_i = 2000;
- dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-
- char *output_filename = "cooling_output.dat";
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- char *output_filename2 = "temperature_output.dat";
- FILE *output_file2 = fopen(output_filename2, "w");
- if (output_file2 == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- char *output_filename3 = "newton_output.dat";
- FILE *output_file3 = fopen(output_filename3, "w");
- if (output_file3 == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- gamma = 5.0 / 3.0;
-
- chemistry_init(params, &us, &internal_const, &chem_data);
- chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- cosmology_init(params, &us, &internal_const, &cosmo);
- cosmology_print(&cosmo);
-
- 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));
-
- cooling_init(params, &us, &internal_const, &cooling);
- cooling_print(&cooling);
-
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- // construct 1d table of cooling rates wrt temperature
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double temp_table[176];
- double element_cooling_table[176];
- double element_print_cooling_table[9 * 176];
- double element_electron_abundance_table[176];
-
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.temperature,
- z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
- cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
- cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
- cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
- construct_1d_print_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
- construct_1d_table_from_3d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
-
-
- //
- // printf("cooling table values \n");
- // for( int j=0; j < 176; j++) {
- // printf(" %.5e", H_plus_He_heat_table[j]+element_cooling_table[j]
- u_ini_cgs = pow(10.0, 14);
- float delta_u;
-
-
- //Compute contributions to cooling rate from different metals
- // int n_t_i = 100;
- // for(int t_i = 0; t_i < n_t_i; t_i++){
- // u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
- // //u = u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
- // u = u_cgs/cooling.internal_energy_scale;
- // 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]);
- // pressure = u*p.rho*(gamma -1.0);
- // p.entropy = pressure/(pow(p.rho,gamma));
- // double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
- // cooling.internal_energy_scale;
-
- // cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
- // H_plus_He_electron_abundance_table,
- // element_print_cooling_table,
- // element_electron_abundance_table,
- // temp_table,
- // z_index, dz, n_h_i, d_n_h, He_i, d_He,
- // &p,&cooling,&cosmo,&internal_const,
- // abundance_ratio);
- // float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
- // &p,&cooling, &cosmo, &internal_const);
- // float temperature_swift = pow(10.0,logT);
-
- // fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- // fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
- // temperature_cgs);
-
- // newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
- // fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
- // //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
- // // \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
- //}
-
- int n_t_i = 100;
- for(int nh_i = 1; nh_i < 7; nh_i++){
- char output_filename4[21];
- sprintf(output_filename4, "%s%d%s", "cooling_output_", nh_i, ".dat");
- FILE *output_file4 = fopen(output_filename4, "w");
- if (output_file4 == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- hydrogen_number_density = pow(10.0,-nh_i)*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo.z,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]);
-
- 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]);
- inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- ratefact = inn_h * (XH / eagle_proton_mass_cgs);
- printf("test cooling hydrogen num dens inn_h %.5e %.5e\n", hydrogen_number_density*cooling.number_density_scale, inn_h);
-
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.temperature,
- z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
- cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
- cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
- cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
- construct_1d_print_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
- construct_1d_table_from_3d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
- n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
- for(int t_i = 0; t_i < n_t_i; t_i++){
- u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
- u = u_cgs/cooling.internal_energy_scale;
- pressure = u*p.rho*(gamma -1.0);
- p.entropy = pressure/(pow(p.rho,gamma));
- double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
- cooling.internal_energy_scale;
-
- double dlambda_du;
- cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u_swift),&dlambda_du,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,NULL,
- abundance_ratio);
- float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
- &p,&cooling, &cosmo, &internal_const);
- float temperature_swift = pow(10.0,logT);
-
- fprintf(output_file4,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
- temperature_cgs);
-
- newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
- fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
- }
- fclose(output_file4);
- }
- fclose(output_file);
- fclose(output_file2);
- fclose(output_file3);
-
- //double dLambdaNet_du, LambdaNet; //, LambdaNext;
- //float x_init, u_eq = 2.0e12;
- //for (int j = 5; j < 6; j++) {
- // float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
- // temp_table, &p, &cooling,
- // &cosmo, &internal_const),
- // x, du;
- // float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
- // LambdaNet = eagle_cooling_rate_1d_table(
- // u_ini, &dLambdaNet_du, H_plus_He_heat_table,
- // H_plus_He_electron_abundance_table, element_cooling_table,
- // element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
- // &internal_const);
- // float u_temp = u_ini + LambdaNet * ratefact * dt;
- // /* RGB removed this **
- // if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
- // &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- // element_cooling_table, element_electron_abundance_table, temp_table, &p,
- // &cooling, &cosmo, &internal_const);
- // if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
- // u_temp = u_ini;
- // } else {
- // u_temp =
- // eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
- // } */
- // if (u_temp > u_eq) {
- // x_init = log(u_temp);
- // } else {
- // x_init = log(u_eq);
- // }
- // x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
- // H_plus_He_electron_abundance_table, element_cooling_table,
- // element_electron_abundance_table, temp_table, &p, &cosmo,
- // &cooling, &internal_const, dt);
- // printf(
- // "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
- // "final %.5e %.5e\n",
- // u_ini, exp(x),
- // eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
- // &cosmo, &internal_const),
- // eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
- // &cosmo, &internal_const));
- //}
-
- free(params);
-
- return 0;
-}
diff --git a/examples/CoolingTest_Alexei/testCooling_script.c b/examples/CoolingTest_Alexei/testCooling_script.c
deleted file mode 100644
index ab4200b67be67e1d45913e56652f9e64913fe3c1..0000000000000000000000000000000000000000
--- a/examples/CoolingTest_Alexei/testCooling_script.c
+++ /dev/null
@@ -1,637 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include <unistd.h>
-#include "cooling.h"
-#include "cooling_struct.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-enum table_index {
- EAGLE_Hydrogen = 0,
- EAGLE_Helium,
- EAGLE_Carbon,
- EAGLE_Nitrogen,
- EAGLE_Oxygen,
- EAGLE_Neon,
- EAGLE_Magnesium,
- EAGLE_Silicon,
- EAGLE_Iron
-};
-
-/*
- * @brief Assign particle density and entropy corresponding to the
- * hydrogen number density and internal energy specified.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param nh Hydrogen number density (cgs units)
- * @param u Internal energy (cgs units)
- */
-void set_quantities(struct part *restrict p,
- const struct unit_system *restrict us,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict internal_const,
- float nh,
- double u){
-
- const float gamma = 5.0/3.0;
- double hydrogen_number_density = nh*pow(units_cgs_conversion_factor(us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo->z,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]);
-
- float scale_factor = 1.0/(1.0+cosmo->z);
- float pressure = (u*pow(scale_factor,3))/cooling->internal_energy_scale*p->rho*(gamma -1.0);
- p->entropy = pressure/(pow(p->rho,gamma));
-}
-
-/*
- * @brief Construct 1d tables from 4d EAGLE tables by
- * interpolating over redshift, hydrogen number density
- * and helium fraction.
- *
- * @param p Particle data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
- * @param temp_table Pointer to 1d interpolated table of temperature values
- * @param H_plus_He_heat_table Pointer to 1d interpolated table of cooling rates
- * due to hydrogen and helium
- * @param H_plus_He_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to hydrogen and helium
- * @param element_electron_abundance_table Pointer to 1d interpolated table
- * of electron abundances due to metals
- * @param element_print_cooling_table Pointer to 1d interpolated table of
- * cooling rates due to each of the metals
- * @param element_cooling_table Pointer to 1d interpolated table of cooling
- * rates due to the contribution of all the metals
- * @param abundance_ratio Pointer to array of ratios of metal abundances to solar
- * @param ub Upper bound in temperature on table construction
- * @param lb Lower bound in temperature on table construction
- */
-void construct_1d_tables_test(const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *restrict internal_const,
- double *temp_table,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_electron_abundance_table,
- double *element_print_cooling_table,
- double *element_cooling_table,
- float *abundance_ratio,
- float *ub, float *lb){
-
- // obtain mass fractions and number density for particle
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, helium fraction, hydrogen number density indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- // construct tables
- //construct_1d_table_from_4d(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.temperature,
- // z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- // cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- //construct_1d_table_from_4d(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- // cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- //construct_1d_table_from_4d(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- // dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
- // cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- //construct_1d_table_from_4d_elements(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- // n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table,abundance_ratio, ub, lb);
- //construct_1d_print_table_from_4d_elements(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- // n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table,abundance_ratio, ub, lb);
- //construct_1d_table_from_3d(
- // p, cooling, cosmo, internal_const,
- // cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
- // n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
-
- if (cosmo->z > cooling->reionisation_redshift) {
- // Photodissociation table
- printf("Eagle testCooling.c photodissociation table redshift reionisation redshift %.5e %.5e\n", cosmo->z, cooling->reionisation_redshift);
- construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.temperature,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,
- He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.H_plus_He_electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
- cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_2d(
- p, cooling, cosmo, internal_const,
- cooling->table.photodissociation_cooling.electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
- printf("Eagle testCooling.c no compton table redshift max redshift %.5e %.5e\n", cosmo->z, cooling->Redshifts[cooling->N_Redshifts - 1]);
- // High redshift table
- construct_1d_table_from_3d(p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.temperature,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.H_plus_He_heating, n_h_i, d_n_h, cooling->N_nH,
- He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.H_plus_He_electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, He_i, d_He, cooling->N_He,
- cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_3d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.metal_heating,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_2d(
- p, cooling, cosmo, internal_const,
- cooling->table.no_compton_cooling.electron_abundance,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- } else {
- // Normal tables
- printf("Eagle testCooling.c normal table redshift %.5e\n", cosmo->z);
- construct_1d_table_from_4d(p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.temperature,
- z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
- cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
- construct_1d_table_from_4d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
- cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
- construct_1d_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table, abundance_ratio, ub, lb);
- construct_1d_print_table_from_4d_elements(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table, abundance_ratio, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, internal_const,
- cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
- n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
- }
-
-}
-
-/*
- * @brief Compare calculation of dlambda/du (gradient of cooling rate,
- * needed for Newton's method) between interpolating 1d tables and
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_dlambda_du(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- // allocate tables
- double H_plus_He_heat_table[cooling->N_Temp];
- double H_plus_He_electron_abundance_table[cooling->N_Temp];
- double temp_table[cooling->N_Temp];
- double element_cooling_table[cooling->N_Temp];
- double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];
- double element_electron_abundance_table[cooling->N_Temp];
- double rate_element_table[cooling->N_Elements+2];
- float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
- double dlambda_du1, dlambda_du2;
-
- // calculate ratio of particle metal abundances to solar
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- // set hydrogen number density and internal energy
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
- // extract hydrogen and helium mass fractions
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, hydrogen number density and helium fraction indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
-
- // construct tables
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling->Temp[0]/eagle_log_10_e;
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-
- // calculate dlambda/du for different values of internal energy
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- p,cooling,cosmo,internal_const,rate_element_table,
- abundance_ratio);
- cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
- &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- p,cooling,cosmo,internal_const,NULL,
- abundance_ratio);
- printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
- }
-}
-
-/*
- * @brief Compare calculating temperature between interpolating 1d tables and
- * 4d tables
- *
- * @param us Units data structure
- * @param p Particle data structure
- * @param xp Extended particle data structure
- * @param internal_const Physical constants data structure
- * @param cooling Cooling function data structure
- * @param cosmo Cosmology data structure
- */
-void compare_temp(
- const struct unit_system *restrict us,
- struct part *restrict p,
- const struct xpart *restrict xp,
- const struct phys_const *restrict internal_const,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo){
-
- // allocate tables
- double H_plus_He_heat_table[cooling->N_Temp];
- double H_plus_He_electron_abundance_table[cooling->N_Temp];
- double temp_table[cooling->N_Temp];
- double element_cooling_table[cooling->N_Temp];
- double element_print_cooling_table[cooling->N_Elements * cooling->N_Temp];
- double element_electron_abundance_table[cooling->N_Temp];
- float *abundance_ratio;
-
- // calculate ratio of particle metal abundances to solar
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- // set hydrogen number density and internal energy
- float nh = 1.0e-1, u = 1.0e9;
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
-
- // extract hydrogen and helium mass fractions
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
-
- // find redshift, hydrogen number density and helium fraction indices and offsets
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
- get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
- get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- // construct tables
- float upper_bound = cooling->Temp[cooling->N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling->Temp[0]/eagle_log_10_e;
- construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
- H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_electron_abundance_table, element_print_cooling_table,
- element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
-
- // calculate temperature for different values of internal energy
- float T1d, T4d, delta_u;
- int nt = 10;
- for(int i = 0; i < nt; i++){
- u = pow(10.0,9 + i);
- set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
- T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- p,cooling, cosmo, internal_const);
- T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, p,cooling, cosmo, internal_const);
- printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
- }
-}
-
-/*
- * @brief Produces contributions to cooling rates for different
- * hydrogen number densities, from different metals,
- * tests 1d and 4d table interpolations produce
- * same results for cooling rate, dlambda/du and temperature.
- *
- * @param pass string "nh" to produce cooling rates for one internal
- * energy, different redshifts
- * @param pass string "metals" to produce contribution to cooling
- * rates from different metals
- * @param pass string "compare_temp" to compare temperature interpolation
- * from 1d and 4d tables
- * @param pass string "compare_dlambda_du" to compare dlambda/du calculation
- * from 1d and 4d tables
- * @param if no string passed, all of the above are performed.
- */
-int main(int argc, char **argv) {
- // Declare relevant structs
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_global_data chem_data;
- struct part p;
- struct xpart xp;
- struct phys_const internal_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "./testCooling.yml";
-
- float nh;
-
- // Read options
- int param, tables = 0;
- float redshift = -1.0, log_10_nh = 100;
- while ((param = getopt(argc, argv, "z:d:t")) != -1)
- switch(param){
- case 'z':
- redshift = atof(optarg);
- break;
- case 'd':
- log_10_nh = atof(optarg);
- break;
- case 't':
- tables = 1;
- break;
- case '?':
- if (optopt == 'z')
- printf ("Option -%c requires an argument.\n", optopt);
- else
- printf ("Unknown option character `\\x%x'.\n",
- optopt);
- error("invalid option(s) to testCooling");
- }
-
- // Read the parameter file
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- // And dump the parameters as used.
- parser_write_params_to_file(params, "used_parameters.yml");
-
- // Init units
- units_init_from_params(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &internal_const);
-
- // Init chemistry
- chemistry_init(params, &us, &internal_const, &chem_data);
- chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
- chemistry_print(&chem_data);
-
- // Init cosmology
- cosmology_init(params, &us, &internal_const, &cosmo);
- cosmology_print(&cosmo);
- if (redshift == -1.0) {
- cosmo.z = 3.0;
- } else {
- cosmo.z = redshift;
- }
-
- // Init cooling
- cooling_init(params, &us, &internal_const, &cooling);
- cooling_print(&cooling);
-
- // Calculate abundance ratios
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
- abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
-
- // Declare 1D tables
- double H_plus_He_heat_table[cooling.N_Temp];
- double H_plus_He_electron_abundance_table[cooling.N_Temp];
- double temp_table[cooling.N_Temp];
- double element_cooling_table[cooling.N_Temp];
- double element_print_cooling_table[cooling.N_Elements * cooling.N_Temp];
- double element_electron_abundance_table[cooling.N_Temp];
-
- // extract mass fractions, calculate table indices and offsets
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
- (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int z_index,He_i,n_h_i;
- float dz,d_He,d_n_h;
- get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
- get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
-
- float upper_bound = cooling.Temp[cooling.N_Temp-1]/eagle_log_10_e;
- float lower_bound = cooling.Temp[0]/eagle_log_10_e;
-
- int nt = 250, n_nh = 6;
- double u = pow(10.0,10);
- //if (argc == 1 || strcmp(argv[1], "nh") == 0){
- // Calculate cooling rates at different densities
- for(int i = 0; i < n_nh; i++){
- // Open files
- char output_filename[21];
- sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density, construct tables
- nh = pow(10.0,-i);
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- temp_table, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_electron_abundance_table,
- element_print_cooling_table,
- element_cooling_table,
- abundance_ratio, &upper_bound, &lower_bound);
- get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
-
- for(int j = 0; j < nt; j++){
- // set internal energy
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
- u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- double dlambda_du;
- float delta_u, cooling_du_dt, logT;
-
- // calculate cooling rates using 1d tables
- if (tables == 1) {
- cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,NULL,
- abundance_ratio);
- logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- &p,&cooling, &cosmo, &internal_const);
- } else {
- // calculate cooling rates using 4d tables
- cooling_du_dt = eagle_metal_cooling_rate(log10(u),
- &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,NULL,
- abundance_ratio);
- logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
- }
- float temperature_swift = pow(10.0,logT);
-
- fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- }
- fclose(output_file);
- }
- //}
- // Calculate contributions from metals to cooling rate
- //if (argc >= 1 || strcmp(argv[1],"metals") == 0) {
- // open file
- char output_filename[21];
- sprintf(output_filename, "%s", "cooling_output.dat");
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density, construct 1d tables
- if (log_10_nh == 100) {
- nh = 1.0e-1;
- } else {
- nh = pow(10.0,log_10_nh);
- }
- //if (argc > 2) nh = pow(10.0,strtod(argv[2],NULL));
- printf("Eagle testcooling.c nh %.5e\n", nh);
- u = pow(10.0,14.0);
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
- construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
- temp_table, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_electron_abundance_table,
- element_print_cooling_table,
- element_cooling_table,
- abundance_ratio, &upper_bound, &lower_bound);
- float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- printf("Eagle testcooling.c inn_h %.5e\n", inn_h);
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
-
- // Loop over internal energy
- for(int j = 0; j < nt; j++){
- set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
- u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
- float delta_u, cooling_du_dt, logT;
-
- // calculate cooling rates using 1d tables
- if (tables == 1) {
- cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_print_cooling_table,
- element_electron_abundance_table,
- temp_table,
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,
- abundance_ratio);
- logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
- &p,&cooling, &cosmo, &internal_const);
- } else {
- // calculate cooling rates using 4d tables
- cooling_du_dt = eagle_print_metal_cooling_rate(
- z_index, dz, n_h_i, d_n_h, He_i, d_He,
- &p,&cooling,&cosmo,&internal_const,
- abundance_ratio);
- logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
- }
- //float temperature_swift = pow(10.0,logT);
- //fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
- fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
- }
- fclose(output_file);
- //}
-
- // compare temperatures and dlambda/du calculated from 1d and 4d tables
- compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
- compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
-
- free(params);
- return 0;
-}
-
-
diff --git a/examples/EAGLE_12/run.sh b/examples/EAGLE_12/run.sh
index 0e35e36d26e110f244265d2e2d3a21381104d655..5674ec36726be8bc1f2d5e63da62835aa211ed7d 100755
--- a/examples/EAGLE_12/run.sh
+++ b/examples/EAGLE_12/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -c -C -t 16 eagle_12.yml 2>&1 | tee output.log
+mpirun -np 4 ../swift_mpi -s -c -C -t 16 -n 100 eagle_12.yml 2>&1 | tee output.log
diff --git a/src/Makefile.am b/src/Makefile.am
index dd6bf25595692b2d0bd2d0a92eb8a6c20909734f..3c7734e554654df747eeb8d1c78ed4c2b85ee130 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -54,7 +54,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
EAGLE_COOLING_SOURCES =
if HAVEEAGLECOOLING
EAGLE_COOLING_SOURCES += cooling/EAGLE/cooling.c \
-cooling/EAGLE/interpolate.c cooling/EAGLE/eagle_cool_tables.c
+cooling/EAGLE/eagle_cool_tables.c
endif
# Common source files
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 6d2b6ce5f647b91c8e778df918341b0539f20046..f4fc179248b45dd79d985518fbfcae4ebeab9d73 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -44,11 +44,12 @@
#include "physical_constants.h"
#include "units.h"
-const float explicit_tolerance = 0.05;
-const float newton_tolerance =
+static const float explicit_tolerance = 0.05;
+static const float newton_tolerance =
1.0e-2; // lower than bisection_tol because using log(u)
-const float bisection_tolerance = 1.0e-6;
-const double bracket_factor = 1.0488088481701; // sqrt(1.1) to match EAGLE
+static const float bisection_tolerance = 1.0e-6;
+static const double bracket_factor = 1.0488088481701; // sqrt(1.1) to match EAGLE
+static float print_cooling_rate_contribution_flag = 0;
/*
* @brief calculates heating due to helium reionization
@@ -73,9 +74,9 @@ eagle_helium_reionization_extraheat(
#endif
/* Helium reionization */
- if (cooling->he_reion == 1) {
+ if (cooling->he_reion_flag == 1) {
double he_reion_erg_pG =
- cooling->he_reion_ev_pH * eagle_ev_to_erg / eagle_proton_mass_cgs;
+ cooling->he_reion_ev_pH / cooling->proton_mass_cgs;
extra_heating += he_reion_erg_pG *
(erf((z - dz - cooling->he_reion_z_center) /
(pow(2., 0.5) * cooling->he_reion_z_sigma)) -
@@ -196,10 +197,10 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
if (z > cooling->Redshifts[cooling->N_Redshifts - 1] ||
z > cooling->reionisation_redshift) {
- T_gam = eagle_cmb_temperature * (1 + z);
+ T_gam = cooling->T_CMB_0 * (1 + z);
- temp_lambda = -eagle_compton_rate *
- (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
+ temp_lambda = -cooling->compton_rate_cgs *
+ (temp - cooling->T_CMB_0 * (1 + z)) * pow((1 + z), 4) *
h_plus_he_electron_abundance / n_h;
cooling_rate += temp_lambda;
if (element_lambda != NULL) element_lambda[1] = temp_lambda;
@@ -284,130 +285,6 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
return cooling_rate;
}
-/*
- * @brief Calculates cooling rate by interpolating precomputed 1d tables
- * of cooling rates which depend only on (log base 10 of) temperature.
- * May be used in bisection scheme when need to perform many table
- * interpolations.
- *
- * @param log_10_u Log base 10 of internal energy
- * @param dlambda_du Pointer to value to be set to rate
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * of change of cooling rate with respect to internal energy
- * @param z_index Redshift index of particle
- * @param dz Redshift offset of particle
- * @param n_h_i Particle hydrogen number density index
- * @param d_n_h Particle hydrogen number density offset
- * @param He_i Particle helium fraction index
- * @param d_He Particle helium fraction offset
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param element_lambda Pointer to array for printing contribution
- * to cooling rate from each of the metals.
- * @param solar_ratio Array of ratios of particle metal abundances
- * to solar metal abundances
- */
-__attribute__((always_inline)) INLINE double eagle_metal_cooling_rate_1d_table(
- double log_10_u, double *dlambda_du, double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table, double *element_cooling_table,
- double *element_electron_abundance_table, double *temp_table,
- 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;
- double cooling_rate = 0.0, temp_lambda;
- float du;
- float h_plus_he_electron_abundance;
-
- int i, temp_i;
- double temp;
- float d_temp;
-
- *dlambda_du = 0.0;
-
- // interpolate to get temperature of particles, find where we are in
- // the temperature table.
- temp = eagle_convert_u_to_temp_1d_table(log_10_u, &du, temp_table, cooling);
- get_index_1d(cooling->Temp, cooling->N_Temp, temp, &temp_i, &d_temp);
-
- /* ------------------ */
- /* Metal-free cooling */
- /* ------------------ */
-
- // contribution to cooling and electron abundance from H, He.
- temp_lambda = interpol_1d_dbl(H_plus_He_heat_table, temp_i, d_temp);
- h_plus_he_electron_abundance =
- interpol_1d_dbl(H_plus_He_electron_abundance_table, temp_i, d_temp);
- cooling_rate += temp_lambda;
- *dlambda_du += (((double)H_plus_He_heat_table[temp_i + 1]) -
- ((double)H_plus_He_heat_table[temp_i])) /
- ((double)du);
- if (element_lambda != NULL) element_lambda[0] = temp_lambda;
-
- /* ------------------ */
- /* Compton cooling */
- /* ------------------ */
-
- // inverse Compton cooling is not in collisional table
- // before reionisation so add now
-
- if (z > cooling->Redshifts[cooling->N_Redshifts - 1] ||
- z > cooling->reionisation_redshift) {
- /* inverse Compton cooling is not in collisional table
- before reionisation so add now */
-
- T_gam = eagle_cmb_temperature * (1 + z);
-
- temp_lambda = -eagle_compton_rate *
- (temp - eagle_cmb_temperature * (1 + z)) * pow((1 + z), 4) *
- h_plus_he_electron_abundance / n_h;
- cooling_rate += temp_lambda;
- if (element_lambda != NULL) element_lambda[1] = temp_lambda;
- }
-
- /* ------------- */
- /* Metal cooling */
- /* ------------- */
-
- // for each element the cooling rate is multiplied by the ratio of H, He
- // electron abundance to solar electron abundance. Note: already multiplied by
- // ratio of particle to solar metal abundances when creating 1d tables.
-
- solar_electron_abundance =
- interpol_1d_dbl(element_electron_abundance_table, temp_i, d_temp);
-
- for (i = 0; i < cooling->N_Elements; i++) {
- temp_lambda = interpol_1d_dbl(element_cooling_table + i * cooling->N_Temp,
- temp_i, d_temp) *
- (h_plus_he_electron_abundance / solar_electron_abundance);
- cooling_rate += temp_lambda;
- *dlambda_du +=
- (((double)element_cooling_table[i * cooling->N_Temp + temp_i + 1]) *
- ((double)H_plus_He_electron_abundance_table[temp_i + 1]) /
- ((double)element_electron_abundance_table[temp_i + 1]) -
- ((double)element_cooling_table[i * cooling->N_Temp + temp_i]) *
- ((double)H_plus_He_electron_abundance_table[temp_i]) /
- ((double)element_electron_abundance_table[temp_i])) /
- ((double)du);
- if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
- }
-
- return cooling_rate;
-}
-
/*
* @brief Wrapper function used to calculate cooling rate and dLambda_du.
* Table indices and offsets for redshift, hydrogen number density and
@@ -441,7 +318,7 @@ __attribute__((always_inline)) INLINE double eagle_cooling_rate(
double lambda_net = 0.0;
lambda_net = eagle_metal_cooling_rate(
- logu / eagle_log_10, dLambdaNet_du, z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ logu / M_LN10, dLambdaNet_du, z_index, dz, n_h_i, d_n_h, He_i, d_He,
p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
if (isnan(lambda_net))
printf(
@@ -452,53 +329,6 @@ __attribute__((always_inline)) INLINE double eagle_cooling_rate(
return lambda_net;
}
-/*
- * @brief Wrapper function used to calculate cooling rate and dLambda_du.
- * Table indices and offsets for redshift, hydrogen number density and
- * helium fraction are passed it so as to compute them only once per particle.
- *
- * @param logu Natural log of internal energy
- * @param dlambda_du Pointer to gradient of cooling rate (set in
- * eagle_metal_cooling_rate)
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * @param z_index Redshift index of particle
- * @param dz Redshift offset of particle
- * @param n_h_i Particle hydrogen number density index
- * @param d_n_h Particle hydrogen number density offset
- * @param He_i Particle helium fraction index
- * @param d_He Particle helium fraction offset
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param abundance_ratio Ratio of element abundance to solar
- */
-__attribute__((always_inline)) INLINE double eagle_cooling_rate_1d_table(
- double logu, double *dLambdaNet_du, double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table, double *element_cooling_table,
- double *element_electron_abundance_table, double *temp_table, int z_index,
- float dz, int n_h_i, float d_n_h, int He_i, float d_He,
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *abundance_ratio) {
- double *element_lambda = NULL, lambda_net = 0.0;
-
- lambda_net = eagle_metal_cooling_rate_1d_table(
- logu / eagle_log_10, dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, p, cooling, cosmo,
- internal_const, element_lambda);
-
- return lambda_net;
-}
-
/*
* @brief Wrapper function used to calculate cooling rate and dLambda_du.
* Writes to file contribution from each element to cooling rate for testing
@@ -530,9 +360,18 @@ __attribute__((always_inline)) INLINE double eagle_print_metal_cooling_rate(
char output_filename[25];
FILE **output_file = malloc((cooling->N_Elements + 2) * sizeof(FILE *));
+ print_cooling_rate_contribution_flag += 1.0/(cooling->N_Elements+2);
for (int element = 0; element < cooling->N_Elements + 2; element++) {
sprintf(output_filename, "%s%d%s", "cooling_element_", element, ".dat");
- output_file[element] = fopen(output_filename, "a");
+ if (print_cooling_rate_contribution_flag < 1) {
+ // If this is the first time we're running this function, overwrite the
+ // output files
+ output_file[element] = fopen(output_filename, "w");
+ print_cooling_rate_contribution_flag += 1.0/(cooling->N_Elements+2);
+ } else {
+ // append to existing files
+ output_file[element] = fopen(output_filename, "a");
+ }
if (output_file == NULL) {
printf("Error opening file!\n");
exit(1);
@@ -552,66 +391,6 @@ __attribute__((always_inline)) INLINE double eagle_print_metal_cooling_rate(
return lambda_net;
}
-/*
- * @brief Wrapper function used to calculate cooling rate and dLambda_du.
- * Prints contribution from each element to cooling rate for testing purposes.
- * Table indices and offsets for redshift, hydrogen number density and
- * helium fraction are passed it so as to compute them only once per particle.
- *
- * @param H_plus_He_heat_table 1D table of heating rates due to H, He
- * @param H_plus_He_electron_abundance_table 1D table of electron abundances due
- * to H,He
- * @param element_cooling_table 1D table of heating rates due to metals
- * @param element_electron_abundance_table 1D table of electron abundances due
- * to metals
- * @param temp_table 1D table of temperatures
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param abundance_ratio Ratio of element abundance to solar
- */
-__attribute__((always_inline)) INLINE double
-eagle_print_metal_cooling_rate_1d_table(
- double *H_plus_He_heat_table, double *H_plus_He_electron_abundance_table,
- double *element_print_cooling_table,
- double *element_electron_abundance_table, double *temp_table,
- 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));
- double u = hydro_get_physical_internal_energy(p, cosmo) *
- cooling->internal_energy_scale;
- double dLambdaNet_du;
-
- char output_filename[25];
- 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_element_", 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;
- lambda_net = eagle_metal_cooling_rate_1d_table(
- log10(u), &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_print_cooling_table,
- element_electron_abundance_table, temp_table, p, cooling, cosmo,
- internal_const, element_lambda);
- for (int j = 0; j < cooling->N_Elements + 2; 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 lambda_net;
-}
-
/*
* @brief Bisection integration scheme used when Newton Raphson method fails to
* converge
@@ -650,7 +429,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
- double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+ double ratefact = inn_h * (XH / cooling->proton_mass_cgs);
// Bracketing
u_lower = u_init;
@@ -760,14 +539,14 @@ __attribute__((always_inline)) INLINE float newton_iter(
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
- double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+ double ratefact = inn_h * (XH / cooling->proton_mass_cgs);
logu_old = logu_init;
logu = logu_old;
int i = 0;
float log_table_bound_high =
- (cooling->Therm[cooling->N_Temp - 1] + 1) / eagle_log_10_e;
- float log_table_bound_low = (cooling->Therm[0] + 1) / eagle_log_10_e;
+ (cooling->Therm[cooling->N_Temp - 1] + 1) / M_LOG10E;
+ float log_table_bound_low = (cooling->Therm[0] + 1) / M_LOG10E;
do /* iterate to convergence */
{
@@ -891,10 +670,10 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
/* 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);
+ ratefact = inn_h * (XH / cooling->proton_mass_cgs);
/* set helium and hydrogen reheating term */
- if (cooling->he_reion == 1)
+ if (cooling->he_reion_flag == 1)
LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
// compute hydrogen number density and helium fraction table indices
@@ -1011,70 +790,6 @@ __attribute__((always_inline)) INLINE float cooling_get_radiated_energy(
return xp->cooling_data.radiated_energy;
}
-/**
- * To do: incorporate into code so that not all tables are read at the beginning
- * @brief Checks the tables that are currently loaded in memory and read
- * new ones if necessary.
- *
- * @param cooling The #cooling_function_data we play with.
- * @param index_z The index along the redshift axis of the tables of the current
- * z.
- */
-inline void eagle_check_cooling_tables(struct cooling_function_data *cooling,
- int index_z) {
-
- /* Do we already have the right table in memory? */
- if (cooling->low_z_index == index_z) return;
-
- if (index_z >= cooling->N_Redshifts) {
-
- error("Missing implementation");
- // Add reading of high-z tables
-
- /* Record the table indices */
- cooling->low_z_index = index_z;
- cooling->high_z_index = index_z;
- } else {
-
- /* Record the table indices */
- cooling->low_z_index = index_z;
- cooling->high_z_index = index_z + 1;
-
- /* Load the damn thing */
- eagle_read_two_tables(cooling);
- }
-}
-
-/**
- * To do: incorporate into code so that not all tables are read at the beginning
- * @brief Common operations performed on the cooling function at a
- * given time-step or redshift.
- *
- * Here we load the cooling tables corresponding to our redshift.
- *
- * @param phys_const The physical constants in internal units.
- * @param us The internal system of units.
- * @param cosmo The current cosmological model.
- * @param cooling The #cooling_function_data used in the run.
- */
-INLINE void cooling_update(const struct phys_const *phys_const,
- const struct unit_system *us,
- const struct cosmology *cosmo,
- struct cooling_function_data *cooling) {
- /* Current redshift */
- const float redshift = cosmo->z;
-
- /* Get index along the redshift index of the tables */
- int index_z;
- float delta_z_table;
- get_redshift_index(redshift, &index_z, &delta_z_table, cooling);
- cooling->index_z = index_z;
- cooling->delta_z_table = delta_z_table;
-
- /* Load the current table (if different from what we have) */
- eagle_check_cooling_tables(cooling, index_z);
-}
-
/**
* @brief Initialises the cooling properties.
*
@@ -1099,7 +814,7 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
parameter_file, "EAGLEChemistry:CalciumOverSilicon");
cooling->sulphur_over_silicon_ratio = parser_get_param_float(
parameter_file, "EAGLEChemistry:SulphurOverSilicon");
- cooling->he_reion =
+ cooling->he_reion_flag =
parser_get_param_float(parameter_file, "EagleCooling:he_reion");
cooling->he_reion_z_center =
parser_get_param_float(parameter_file, "EagleCooling:he_reion_z_center");
@@ -1108,15 +823,16 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
cooling->he_reion_ev_pH =
parser_get_param_float(parameter_file, "EagleCooling:he_reion_ev_pH");
+ /* convert to cgs */
+ cooling->he_reion_ev_pH *= phys_const->const_electron_volt *
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY);
+
// read in cooling tables
GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname, cooling);
cooling->table = eagle_readtable(cooling->cooling_table_path, cooling);
- // allocate array for solar abundances
- cooling->solar_abundances = malloc(cooling->N_Elements * sizeof(double));
-
// compute conversion factors
cooling->internal_energy_scale =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) /
@@ -1128,6 +844,31 @@ INLINE void cooling_init_backend(struct swift_params *parameter_file,
units_cgs_conversion_factor(us, UNIT_CONV_MASS);
cooling->temperature_scale =
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+ cooling->proton_mass_cgs = phys_const->const_proton_mass *
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ cooling->T_CMB_0 = phys_const->const_T_CMB_0 *
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
+
+ /* Compute the coefficient at the front of the Compton cooling expression */
+ const double radiation_constant =
+ 4. * phys_const->const_stefan_boltzmann / phys_const->const_speed_light_c;
+ const double compton_coefficient =
+ 4. * radiation_constant * phys_const->const_thomson_cross_section *
+ phys_const->const_boltzmann_k /
+ (phys_const->const_electron_mass * phys_const->const_speed_light_c);
+ const float dimension_coefficient[5] = {1, 2, -3, 0, -5};
+
+ /* This should be ~1.0178085e-37 g cm^2 s^-3 K^-5 */
+ const double compton_coefficient_cgs =
+ compton_coefficient *
+ units_general_cgs_conversion_factor(us, dimension_coefficient);
+
+ /* And now the Compton rate */
+ cooling->compton_rate_cgs =
+ compton_coefficient_cgs * cooling->T_CMB_0 *
+ cooling->T_CMB_0 * cooling->T_CMB_0 * cooling->T_CMB_0;
+
}
/**
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index f1ad66b055861cd80c5ab9b1d99799449ac5ed79..2eb887ae51377e6922699d16c352fce11d6392aa 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -32,15 +32,6 @@
#include "physical_constants.h"
#include "units.h"
-enum hdf5_allowed_types {
- hdf5_short,
- hdf5_int,
- hdf5_long,
- hdf5_float,
- hdf5_double,
- hdf5_char
-};
-
double eagle_helium_reionization_extraheat(
double, double, const struct cooling_function_data *restrict);
@@ -50,35 +41,17 @@ double eagle_metal_cooling_rate(double, double *, int, float, int, float, int,
const struct cosmology *restrict,
const struct phys_const *, double *, float *);
-double eagle_metal_cooling_rate_1d_table(
- double, double *, double *, double *, double *, double *, double *,
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *, double *);
-
double eagle_cooling_rate(double, double *, int, float, int, float, int, float,
const struct part *restrict,
const struct cooling_function_data *restrict,
const struct cosmology *restrict,
const struct phys_const *, float *);
-double eagle_cooling_rate_1d_table(double, double *, double *, double *,
- double *, double *, double *, int, float,
- int, float, int, float,
- const struct part *restrict,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *, float *);
-
double eagle_print_metal_cooling_rate(
int, float, int, float, int, float, const struct part *restrict,
const struct cooling_function_data *restrict,
const struct cosmology *restrict, const struct phys_const *, float *);
-double eagle_print_metal_cooling_rate_1d_table(
- double *, double *, double *, double *, double *,
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *);
-
float bisection_iter(float, double, int, float, int, float, int, float, float,
struct part *restrict, const struct cosmology *restrict,
const struct cooling_function_data *restrict,
@@ -116,11 +89,6 @@ void cooling_first_init_part(const struct phys_const *restrict,
float cooling_get_radiated_energy(const struct xpart *restrict);
-void eagle_check_cooling_tables(struct cooling_function_data *, int);
-
-void cooling_update(const struct phys_const *, const struct unit_system *,
- const struct cosmology *, struct cooling_function_data *);
-
void cooling_init_backend(struct swift_params *, const struct unit_system *,
const struct phys_const *,
struct cooling_function_data *);
diff --git a/src/cooling/EAGLE/cooling_read_table.h b/src/cooling/EAGLE/cooling_read_table.h
deleted file mode 100644
index 82bf92bdd54020bce12c9ffc8a42b19b0ea400fa..0000000000000000000000000000000000000000
--- a/src/cooling/EAGLE/cooling_read_table.h
+++ /dev/null
@@ -1,132 +0,0 @@
-#ifndef SWIFT_COOLING_EAGLE_READ_TABLES_H
-#define SWIFT_COOLING_EAGLE_READ_TABLES_H
-
-/*
- * @brief Reads cooling header data into cooling data structure
- *
- * @param fname Filepath
- * @param cooling Cooling data structure
- */
-INLINE void ReadCoolingHeader(char *fname,
- struct cooling_function_data *cooling) {
- 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) {
- error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
- }
-
- dataset =
- H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
- 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_density_bins", H5P_DEFAULT);
- 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_helium_fractions", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_He);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_SolarAbundances);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_Elements);
- status = H5Dclose(dataset);
-
- /* allocate arrays for cooling table header */
- // allocate_header_arrays();
-
- cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
- cooling->nH = malloc(cooling->N_Temp * sizeof(float));
- cooling->HeFrac = malloc(cooling->N_Temp * sizeof(float));
- cooling->SolarAbundances = malloc(cooling->N_Temp * sizeof(float));
-
- /* fill the arrays */
- dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
- 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", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->nH);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
- H5P_DEFAULT);
- 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",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->SolarAbundances);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->Therm);
- status = H5Dclose(dataset);
-
- char element_names[cooling->N_Elements][element_name_length];
- hsize_t string_length = element_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", H5P_DEFAULT);
- 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",
- // H5P_DEFAULT);
- // 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);
-}
-
-#endif /* SWIFT_COOLING_EAGLE_READ_TABLES_H */
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index da157ad8db558a162bc70c753d59c7e91e9b303c..8b86239e03cb91621d25348067d3ffbb6fa700e7 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -18,31 +18,11 @@
******************************************************************************/
#ifndef SWIFT_COOLING_STRUCT_EAGLE_H
#define SWIFT_COOLING_STRUCT_EAGLE_H
-#include <stdbool.h>
-#include <time.h>
// EAGLE defined constants
#define eagle_element_name_length 20
-#define eagle_cmb_temperature 2.728
-#define eagle_compton_rate 1.0178e-37 * 2.728 * 2.728 * 2.728 * 2.728
#define eagle_metal_cooling_on 1
#define eagle_max_iterations 15
-#define eagle_proton_mass_cgs 1.6726e-24
-#define eagle_ev_to_erg 1.60217646e-12
-#define eagle_log_10 2.30258509299404
-#define eagle_log_10_e 0.43429448190325182765
-
-/*
- * @brief struct containing radiative and heating rates
- */
-// struct radiative_rates{
-// float Cooling[NUMBER_OF_CTYPES];
-// float Heating[NUMBER_OF_HTYPES];
-// float TotalCoolingRate;
-// float TotalHeatingRate;
-// float CoolingTime; /*Cooling time in Myr*/
-// float HeatingTime; /*Heating time in Myr*/
-//};
/*
* @brief struct containing cooling tables independent of redshift
@@ -99,50 +79,59 @@ struct cooling_function_data {
/* Cooling table variables */
struct eagle_cooling_table table;
+ /* Size of table dimensions */
int N_Redshifts;
int N_nH;
int N_Temp;
int N_He;
+
+ /* Number of metals and solar abundances tracked in EAGLE */
int N_Elements;
int N_SolarAbundances;
+ /* Arrays of grid values in tables */
float *Redshifts;
float *nH;
float *Temp;
float *HeFrac;
float *Therm;
+
+ /* Array of values of solar metal and electron abundance */
float *SolarAbundances;
float *SolarElectronAbundance;
- float *ElementAbundance_SOLARM1;
- double *solar_abundances;
+
+
+ /* Arrays containing names of tracked elements. Used for
+ * reading in the contributions to the cooling rate from
+ * each element */
char **ElementNames;
char **SolarAbundanceNames;
- int *ElementNamePointers;
- int *SolarAbundanceNamePointers;
-
- /*! Constant multiplication factor for time-step criterion */
- float cooling_tstep_mult;
- float Redshift;
- float min_energy;
- float cooling_rate;
+
+ /* Normalisation constants that are frequently used */
double internal_energy_scale;
double number_density_scale;
double temperature_scale;
double power_scale;
+
+ /* filepath to EAGLE cooling tables */
char cooling_table_path[500];
+
+ /* Some constants read in from yml file relevant to EAGLE cooling */
float reionisation_redshift;
float calcium_over_silicon_ratio;
float sulphur_over_silicon_ratio;
- int he_reion;
+ /* Helium reionisation parameters */
+ int he_reion_flag;
float he_reion_ev_pH;
float he_reion_z_center;
float he_reion_z_sigma;
- int index_z, low_z_index, high_z_index;
- float delta_z_table;
+ /* Proton mass in cgs */
+ double proton_mass_cgs;
+ double T_CMB_0;
+ double compton_rate_cgs;
- double delta_u;
};
/**
diff --git a/src/cooling/EAGLE/eagle_cool_tables.c b/src/cooling/EAGLE/eagle_cool_tables.c
index 9b4bd6b17ba43c2112368afd30e35ff16dfee76f..4aa8d8496095ebd703d4ecf5057c7037223aca25 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.c
+++ b/src/cooling/EAGLE/eagle_cool_tables.c
@@ -36,19 +36,6 @@
#include "error.h"
#include "interpolate.h"
-/*
- * @brief String assignment function from EAGLE
- *
- * @param s String
- */
-char *mystrdup(const char *s) {
- char *p;
-
- p = (char *)malloc((strlen(s) + 1) * sizeof(char));
- strcpy(p, s);
- return p;
-}
-
/*
* @brief Reads in EAGLE table of redshift values
*
@@ -134,10 +121,14 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
cooling->HeFrac = malloc(cooling->N_He * sizeof(float));
cooling->SolarAbundances = malloc(cooling->N_SolarAbundances * sizeof(float));
cooling->Therm = malloc(cooling->N_Temp * sizeof(float));
- cooling->ElementNames =
- malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
- cooling->SolarAbundanceNames = malloc(
- cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
+
+ cooling->ElementNames = malloc(cooling->N_Elements * sizeof(char *));
+ for (i = 0; i < cooling->N_Elements; i++)
+ cooling->ElementNames[i] = malloc(eagle_element_name_length * sizeof(char));
+
+ cooling->SolarAbundanceNames = malloc(cooling->N_SolarAbundances * sizeof(char *));
+ for (i = 0; i < cooling->N_SolarAbundances; i++)
+ cooling->SolarAbundanceNames[i] = malloc(eagle_element_name_length * sizeof(char));
// read in values for each of the arrays
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
@@ -181,7 +172,8 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
H5Tclose(datatype);
for (i = 0; i < cooling->N_Elements; i++)
- cooling->ElementNames[i] = mystrdup(element_names[i]);
+ strcpy(cooling->ElementNames[i],element_names[i]);
+
char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
@@ -196,7 +188,7 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
H5Tclose(datatype);
for (i = 0; i < cooling->N_SolarAbundances; i++)
- cooling->SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
+ strcpy(cooling->SolarAbundanceNames[i],element_names[i]);
status = H5Fclose(tempfile_id);
@@ -355,7 +347,7 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- printf("eagle_cool_tables.h done reading in no compton table\n");
+ message("done reading in no compton table");
return cooling_table;
}
@@ -490,7 +482,7 @@ struct cooling_tables_redshift_invariant get_collisional_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- printf("eagle_cool_tables.h done reading in collisional table\n");
+ message("done reading in collisional table");
return cooling_table;
}
@@ -639,7 +631,7 @@ struct cooling_tables_redshift_invariant get_photodis_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- printf("eagle_cool_tables.h done reading in photodissociation table\n");
+ message("done reading in photodissociation table");
return cooling_table;
}
@@ -793,7 +785,7 @@ struct cooling_tables get_cooling_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- printf("eagle_cool_tables.h done reading in general cooling table\n");
+ message("done reading in general cooling table");
return cooling_table;
}
@@ -820,179 +812,4 @@ struct eagle_cooling_table eagle_readtable(
return table;
}
-/*
- * @brief Work in progress: read in only two tables for given redshift
- *
- * @param cooling_table_path Directory containing cooling tables
- * @param cooling Cooling function data structure
- * @param filename1 filename2 Names of the two tables we need to read
- */
-struct cooling_tables get_two_cooling_tables(
- char *cooling_table_path,
- const struct cooling_function_data *restrict cooling, char *filename1,
- char *filename2) {
-
- struct cooling_tables 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;
- float *electron_abundance;
- float *temperature;
- float *he_net_cooling_rate;
- float *he_electron_abundance;
-
- // allocate temporary arrays (needed to change order of dimensions
- // of arrays)
- net_cooling_rate =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- // allocate arrays that the values will be assigned to
- cooling_table.metal_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
- cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(
- cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.temperature =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_electron_abundance =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- // repeat for both tables we need to read
- for (int file = 0; file < 2; file++) {
- if (file == 0) {
- sprintf(fname, "%s%s.hdf5", cooling_table_path, filename1);
- } else {
- sprintf(fname, "%s%s.hdf5", cooling_table_path, filename2);
- }
- file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
- if (file_id < 0) {
- error("[GetCoolingTables()]: unable to open file %s\n", fname);
- }
-
- // read in cooling rates due to metals
- 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 (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(
- file, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
- cooling->N_Temp, cooling->N_Elements);
- cooling_table.metal_heating[cooling_index] =
- -net_cooling_rate[table_index];
- }
- }
- }
-
- // read in cooling rates due to hydrogen and helium, H + He electron
- // abundances, temperatures
- 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++) {
- 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(file, k, i, j, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, 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];
- cooling_table.temperature[cooling_index] =
- log10(temperature[table_index]);
- }
- }
- }
-
- // read in electron densities due to metals
- 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++) {
- 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(file, j, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
- cooling_table.electron_abundance[cooling_index] =
- electron_abundance[table_index];
- }
- }
-
- 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 general cooling table\n");
-
- return cooling_table;
-}
-
-/*
- * @brief Work in progress: constructs the data structure containting cooling
- * tables bracketing the redshift of a particle.
- *
- * @param cooling_table_path Directory containing cooling table
- * @param cooling Cooling data structure
- */
-void eagle_read_two_tables(struct cooling_function_data *restrict cooling) {
-
- struct eagle_cooling_table table;
-
- table.element_cooling =
- get_cooling_table(cooling->cooling_table_path, cooling);
- cooling->table = table;
-}
-
#endif
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index 760ed59487f549e026cd6a63ae533a785ba868b1..3899d5396b49efad5a1361bab41c2d156ba0d7e5 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -1,1020 +1,57 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
#ifndef SWIFT_EAGLE_COOL_TABLES_H
#define SWIFT_EAGLE_COOL_TABLES_H
+/**
+ * @file src/cooling/EAGLE/cooling.h
+ * @brief EAGLE cooling function
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
#include <hdf5.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "cooling_struct.h"
+#include "interpolate.h"
#include "error.h"
-/*
- * @brief String assignment function from EAGLE
- *
- * @param s String
- */
-inline char *mystrdup(const char *s) {
- char *p;
-
- p = (char *)malloc((strlen(s) + 1) * sizeof(char));
- strcpy(p, s);
- return p;
-}
-
-int row_major_index_2d(int, int, int, int);
-int row_major_index_3d(int, int, int, int, int, int);
-int row_major_index_4d(int, int, int, int, int, int, int, int);
-
-/*
- * @brief Reads in EAGLE table redshift values
- *
- * @param cooling Cooling data structure
- */
-inline void GetCoolingRedshifts(struct cooling_function_data *cooling) {
- FILE *infile;
-
- int i = 0;
-
- char buffer[500], redfilename[500];
-
- sprintf(redfilename, "%s/redshifts.dat", cooling->cooling_table_path);
- 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 *)malloc(cooling->N_Redshifts * sizeof(float));
-
- while (fscanf(infile, "%s", buffer) != EOF) {
- cooling->Redshifts[i] = atof(buffer);
- i += 1;
- }
- }
- fclose(infile);
-}
-
-/*
- * @brief Reads in EAGLE cooling table header
- *
- * @param fname Filepath
- * @param cooling Cooling data structure
- */
-inline void ReadCoolingHeader(char *fname,
- struct cooling_function_data *cooling) {
- 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) {
- error("[ReadCoolingHeader()]: unable to open file %s\n", fname);
- }
-
- dataset =
- H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins", H5P_DEFAULT);
- 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_density_bins", H5P_DEFAULT);
- 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_helium_fractions", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_He);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_SolarAbundances);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- &cooling->N_Elements);
- status = H5Dclose(dataset);
-
- /* allocate arrays for cooling table header */
- // allocate_header_arrays();
-
- cooling->Temp = malloc(cooling->N_Temp * sizeof(float));
- cooling->nH = malloc(cooling->N_nH * sizeof(float));
- cooling->HeFrac = malloc(cooling->N_He * sizeof(float));
- cooling->SolarAbundances = malloc(cooling->N_SolarAbundances * sizeof(float));
- cooling->Therm = malloc(cooling->N_Temp * sizeof(float));
- cooling->ElementNames =
- malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
- cooling->SolarAbundanceNames = malloc(
- cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
-
- /* fill the arrays */
- dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
- 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", H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->nH);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins",
- H5P_DEFAULT);
- 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",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->SolarAbundances);
- status = H5Dclose(dataset);
-
- dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins",
- H5P_DEFAULT);
- status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
- cooling->Therm);
- status = H5Dclose(dataset);
-
- char element_names[cooling->N_Elements][eagle_element_name_length];
- hsize_t string_length = eagle_element_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", H5P_DEFAULT);
- status =
- H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
- status = H5Dclose(dataset);
- H5Tclose(datatype);
-
- for (i = 0; i < cooling->N_Elements; i++)
- cooling->ElementNames[i] = mystrdup(element_names[i]);
-
- char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
-
- /* assumed solar abundances used in constructing the tables, and corresponding
- * names */
- datatype = H5Tcopy(H5T_C_S1);
- status = H5Tset_size(datatype, string_length);
- dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names", H5P_DEFAULT);
- 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);
-}
-
-/*
- * @brief Get the cooling table for photoionized cooling (before redshift ~9)
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables_redshift_invariant get_no_compt_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;
- float *electron_abundance;
- float *temperature;
- float *he_net_cooling_rate;
- float *he_electron_abundance;
-
- net_cooling_rate =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(
- cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(
- cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_electron_abundance = (float *)malloc(
- cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
-
- sprintf(fname, "%sz_8.989nocompton.hdf5", cooling_table_path);
-
- file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
- // printf("GetNoCompTable Redshift 1 %ld %s\n", (long int)file_id, fname);
- // fflush(stdout);
-
- /* 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++) {
- 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, k, j, specs, 1, cooling->N_nH, cooling->N_Temp,
- cooling->N_Elements); // Redshift invariant table!!!
- cooling_table.metal_heating[cooling_index] =
- -net_cooling_rate[table_index];
- }
- }
- }
-
- /* 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++) {
- 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, k, i, j, 1, cooling->N_nH, 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++) {
- 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_table.electron_abundance[cooling_index] =
- electron_abundance[table_index];
- }
- }
-
- status = H5Fclose(file_id);
-
- // cooling_table.metals_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 no compton table\n");
-
- return cooling_table;
-}
-
-/*
- * @brief Get the cooling table for collisional cooling (before reionisation)
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables_redshift_invariant get_collisional_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, 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;
-}
-
-/*
- * @brief Get the cooling table for photodissociation
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-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;
- float *electron_abundance;
- float *temperature;
- float *he_net_cooling_rate;
- float *he_electron_abundance;
-
- net_cooling_rate =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- cooling_table.metal_heating = (float *)malloc(
- cooling->N_Elements * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_heating = (float *)malloc(
- cooling->N_He * cooling->N_Temp * cooling->N_nH * sizeof(float));
- 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);
-
- 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++) {
- 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_3d(
- k, j, specs, cooling->N_nH, cooling->N_Temp,
- cooling->N_Elements); // Redshift invariant table!!!
- cooling_table.metal_heating[cooling_index] =
- -net_cooling_rate[table_index];
- }
- }
- }
-
- /* 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++) {
- 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_3d(k, i, j, cooling->N_nH, 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++) {
- 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_2d(
- j, i, cooling->N_nH, cooling->N_Temp); // Redshift invariant table!!!
- cooling_table.electron_abundance[cooling_index] =
- electron_abundance[table_index];
- }
- }
-
- 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 photodissociation table\n");
- return cooling_table;
-}
-
-inline struct cooling_tables get_two_cooling_tables(
- char *cooling_table_path,
- const struct cooling_function_data *restrict cooling, char *filename1,
- char *filename2) {
-
- struct cooling_tables 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;
- float *electron_abundance;
- float *temperature;
- float *he_net_cooling_rate;
- float *he_electron_abundance;
-
- net_cooling_rate =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- cooling_table.metal_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
- cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(
- cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.temperature =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_electron_abundance =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- /* For normal elements */
- for (int file = 0; file < 2; file++) {
- if (file == 0) {
- sprintf(fname, "%s%s.hdf5", cooling_table_path, filename1);
- } else {
- sprintf(fname, "%s%s.hdf5", cooling_table_path, filename2);
- }
- file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
- if (file_id < 0) {
- error("[GetCoolingTables()]: unable to open file %s\n", fname);
- }
-
- 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 (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(
- file, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
- cooling->N_Temp, cooling->N_Elements);
- cooling_table.metal_heating[cooling_index] =
- -net_cooling_rate[table_index];
- }
- }
- }
-
- /* 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++) {
- 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(file, k, i, j, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, 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];
- 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++) {
- 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(file, j, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
- cooling_table.electron_abundance[cooling_index] =
- electron_abundance[table_index];
- }
- }
-
- 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 general cooling table\n");
-
- return cooling_table;
-}
-
-/*
- * @brief Get the cooling tables dependent on redshift
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-
-inline struct cooling_tables get_cooling_table(
- char *cooling_table_path,
- const struct cooling_function_data *restrict cooling) {
-
- struct cooling_tables 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;
- float *electron_abundance;
- float *temperature;
- float *he_net_cooling_rate;
- float *he_electron_abundance;
-
- net_cooling_rate =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- electron_abundance =
- (float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
- temperature = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_net_cooling_rate = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- he_electron_abundance = (float *)malloc(cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- cooling_table.metal_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_Elements *
- cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.electron_abundance = (float *)malloc(
- cooling->N_Redshifts * cooling->N_Temp * cooling->N_nH * sizeof(float));
- cooling_table.temperature =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_heating =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
- cooling_table.H_plus_He_electron_abundance =
- (float *)malloc(cooling->N_Redshifts * cooling->N_He * cooling->N_Temp *
- cooling->N_nH * sizeof(float));
-
- /* For normal elements */
- for (int z_index = 0; z_index < cooling->N_Redshifts; z_index++) {
- sprintf(fname, "%sz_%1.3f.hdf5", cooling_table_path,
- cooling->Redshifts[z_index]);
- file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
-
- if (file_id < 0) {
- error("[GetCoolingTables()]: unable to open file %s\n", fname);
- }
-
- 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 (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, i, j, specs, cooling->N_Redshifts, cooling->N_nH,
- cooling->N_Temp, cooling->N_Elements);
- cooling_table.metal_heating[cooling_index] =
- -net_cooling_rate[table_index];
- }
- }
- }
-
- /* 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++) {
- 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, k, i, j, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, 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];
- 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++) {
- 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_table.electron_abundance[cooling_index] =
- electron_abundance[table_index];
- }
- }
-
- 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 general cooling table\n");
-
- return cooling_table;
-}
-
-/*
- * @brief Constructs the data structure containting all the cooling tables
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-inline void eagle_read_two_tables(
- struct cooling_function_data *restrict cooling) {
-
- struct eagle_cooling_table table;
-
- table.element_cooling =
- get_cooling_table(cooling->cooling_table_path, cooling);
- cooling->table = table;
-}
-
-/*
- * @brief Constructs the data structure containting all the cooling tables
- *
- * @param cooling_table_path Filepath
- * @param cooling Cooling data structure
- */
-inline struct eagle_cooling_table eagle_readtable(
- char *cooling_table_path,
- const struct cooling_function_data *restrict cooling) {
-
- struct eagle_cooling_table table;
-
- table.no_compton_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);
-
- return table;
-}
-
-/*
- * @brief Finds the element index for the corresponding element in swift
- *
- * @param element_name Element string we want to match to element index
- * @param cooling Cooling data structure
- */
-inline int element_index(char *element_name,
- const struct cooling_function_data *restrict cooling) {
- int i;
-
- for (i = 0; i < cooling->N_Elements; i++)
- if (strcmp(cooling->ElementNames[i], element_name) == 0) return i;
-
- /* element not found */
- return -1;
-}
-
-/*
- * @brief Finds the element index for the corresponding element in EAGLE
- *
- * @param element_name Element string we want to match to element index
- * @param size Number of elements tracked in EAGLE
- * @param cooling Cooling data structure
- */
-inline int get_element_index(char *table[20], int size, char *element_name) {
- int i;
-
- for (i = 0; i < size; i++)
- if (strcmp(table[i], element_name) == 0) return i;
-
- /* element not found */
- return -1;
-}
-
-/*
- * @brief Makes an array of element names which are tracked in the EAGLE tables
- *
- * @param cooling Cooling data structure
- */
-inline void MakeNamePointers(struct cooling_function_data *cooling) {
- int i, j, sili_index = 0;
- char ElementNames[cooling->N_Elements][eagle_element_name_length];
+void GetCoolingRedshifts(struct cooling_function_data *);
- /* This is ridiculous, way too many element name arrays. Needs to be changed
- */
- // ElementNames =
- // malloc(cooling->N_Elements*eagle_element_name_length*sizeof(char));
- strcpy(ElementNames[0], "Hydrogen");
- strcpy(ElementNames[1], "Helium");
- strcpy(ElementNames[2], "Carbon");
- strcpy(ElementNames[3], "Nitrogen");
- strcpy(ElementNames[4], "Oxygen");
- strcpy(ElementNames[5], "Neon");
- strcpy(ElementNames[6], "Magnesium");
- strcpy(ElementNames[7], "Silicon");
- strcpy(ElementNames[8], "Iron");
+void ReadCoolingHeader(char *, struct cooling_function_data *);
- cooling->ElementNamePointers = malloc(cooling->N_Elements * sizeof(int));
- cooling->SolarAbundanceNamePointers =
- malloc(cooling->N_Elements * sizeof(int));
+struct cooling_tables_redshift_invariant get_no_compt_table(
+ char *, const struct cooling_function_data *restrict);
- for (i = 0; i < cooling->N_Elements; i++) {
- if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
- }
+struct cooling_tables_redshift_invariant get_collisional_table(
+ char *, const struct cooling_function_data *restrict);
- for (i = 0; i < cooling->N_Elements; i++) {
- cooling->SolarAbundanceNamePointers[i] = -999;
- cooling->ElementNamePointers[i] = -999;
+struct cooling_tables_redshift_invariant get_photodis_table(
+ char *, const struct cooling_function_data *restrict);
- for (j = 0; j < cooling->N_SolarAbundances; j++) {
- if (strcmp(cooling->ElementNames[i], cooling->SolarAbundanceNames[j]) ==
- 0)
- cooling->SolarAbundanceNamePointers[i] = j;
- }
+struct cooling_tables get_cooling_table(
+ char *, const struct cooling_function_data *restrict);
- if (strcmp(cooling->ElementNames[i], "Sulphur") == 0 ||
- strcmp(cooling->ElementNames[i], "Calcium") ==
- 0) /* These elements are tracked! */
- cooling->ElementNamePointers[i] = -1 * sili_index;
- else {
- for (j = 0; j < cooling->N_Elements; j++) {
- if (strcmp(cooling->ElementNames[i], ElementNames[j]) == 0)
- cooling->ElementNamePointers[i] = j;
- }
- }
- }
-}
+struct eagle_cooling_table eagle_readtable(
+ char *, const struct cooling_function_data *restrict);
#endif
diff --git a/src/cooling/EAGLE/interpolate.h b/src/cooling/EAGLE/interpolate.h
index 151bc735067d55481befa5ee5f62815b93174c27..ed82ee26e5a02e9663b9dfe1524b506a1871eb75 100644
--- a/src/cooling/EAGLE/interpolate.h
+++ b/src/cooling/EAGLE/interpolate.h
@@ -21,7 +21,7 @@
/**
* @file src/cooling/EAGLE/interpolate.h
- * @brief EAGLE interpolation function declarations
+ * @brief Interpolation functions for EAGLE tables
*/
/* Config parameters. */
@@ -36,6 +36,7 @@
/* Local includes. */
#include "chemistry.h"
#include "cooling_struct.h"
+#include "eagle_cool_tables.h"
#include "error.h"
#include "hydro.h"
#include "io_properties.h"
@@ -44,90 +45,415 @@
#include "physical_constants.h"
#include "units.h"
-int row_major_index_2d(int, int, int, int);
+static int get_redshift_index_first_call = 0;
+static int get_redshift_index_previous = -1;
-int row_major_index_3d(int, int, int, int, int, int);
+static const float EPS = 1.e-4;
-int row_major_index_4d(int, int, int, int, int, int, int, int);
+/**
+ * @brief Returns the 1d index of element with 2d indices i,j
+ * from a flattened 2d array in row major order
+ *
+ * @param i, j Indices of element of interest
+ * @param nx, ny Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
+ int nx, int ny) {
+ return i * ny + j;
+}
-void get_index_1d(float *, int, double, int *, float *);
+/**
+ * @brief Returns the 1d index of element with 3d indices i,j,k
+ * from a flattened 3d array in row major order
+ *
+ * @param i, j, k Indices of element of interest
+ * @param nx, ny, nz Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
+ int k, int nx,
+ int ny, int nz) {
+ return i * ny * nz + j * nz + k;
+}
-void get_redshift_index(float, int *, float *,
- const struct cooling_function_data *restrict);
+/**
+ * @brief Returns the 1d index of element with 4d indices i,j,k,l
+ * from a flattened 4d array in row major order
+ *
+ * @param i, j, k, l Indices of element of interest
+ * @param nx, ny, nz, nw Sizes of array dimensions
+ */
+__attribute__((always_inline)) INLINE int row_major_index_4d(int i, int j,
+ int k, int l,
+ int nx, int ny,
+ int nz, int nw) {
+ return i * ny * nz * nw + j * nz * nw + k * nw + l;
+}
-float interpol_1d(float *, int, float);
+/*
+ * @brief 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.
+ *
+ * @param table Pointer to table of values
+ * @param ntable Size of the table
+ * @param x Value whose position within the array we are interested in
+ * @param i Pointer to the index whose corresponding table value
+ * is the greatest value in the table less than x
+ * @param dx Pointer to offset of x within table cell
+ */
+__attribute__((always_inline)) INLINE void get_index_1d(float *table,
+ int ntable, double x,
+ int *i, float *dx) {
-double interpol_1d_dbl(double *, int, float);
+ float dxm1 = (float)(ntable - 1) / (table[ntable - 1] - table[0]);
-float interpol_2d(float *, int, int, float, float, int, int, double *,
- double *);
+ 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;
+ }
+}
-double interpol_2d_dbl(double *, int, int, double, double, int, int, double *,
- double *);
+/*
+ * @brief Returns the position i of a value z in the redshift table
+ * and computes the displacement dz needed for the interpolation.
+ * Since the redshift table is not evenly spaced, compare z with each
+ * table value in decreasing order starting with the previous redshift index
+ *
+ * @param z Redshift whose position within the redshift array we are interested
+ * in
+ * @param z_index Pointer to the index whose corresponding redshift
+ * is the greatest value in the redshift table less than x
+ * @param dz Pointer to offset of z within redshift cell
+ * @param cooling Pointer to cooling structure containing redshift table
+ */
+__attribute__((always_inline)) INLINE void get_redshift_index(
+ float z, int *z_index, float *dz,
+ const struct cooling_function_data *restrict cooling) {
+ int i, iz;
-float interpol_3d(float *, int, int, int, float, float, float, int, int, int,
- double *, double *);
+ if (get_redshift_index_first_call == 0) {
+ get_redshift_index_first_call = 1;
+ get_redshift_index_previous = cooling->N_Redshifts - 2;
-float interpol_4d(float *, int, int, int, int, float, float, float, float, int,
- int, int, int, double *, double *);
+ /* this routine assumes cooling_redshifts table is in increasing order. Test
+ * this. */
+ for (i = 0; i < cooling->N_Redshifts - 2; i++)
+ if (cooling->Redshifts[i + 1] < cooling->Redshifts[i]) {
+ error("[get_redshift_index]: table should be in increasing order\n");
+ }
+ }
-void construct_1d_table_from_2d(const struct part *restrict,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *, float *, int, float,
- int, int, double *, float *, float *);
+ /* before the earliest redshift or before hydrogen reionization, flag for
+ * collisional cooling */
+ if (z > cooling->reionisation_redshift) {
+ *z_index = cooling->N_Redshifts;
+ *dz = 0.0;
+ }
+ /* from reionization use the cooling tables */
+ else if (z > cooling->Redshifts[cooling->N_Redshifts - 1] &&
+ z <= cooling->reionisation_redshift) {
+ *z_index = cooling->N_Redshifts + 1;
+ *dz = 0.0;
+ }
+ /* at the end, just use the last value */
+ else if (z <= cooling->Redshifts[0]) {
+ *z_index = 0;
+ *dz = 0.0;
+ } else {
+ /* start at the previous index and search */
+ for (iz = get_redshift_index_previous; iz >= 0; iz--) {
+ if (z > cooling->Redshifts[iz]) {
+ *dz = (z - cooling->Redshifts[iz]) /
+ (cooling->Redshifts[iz + 1] - cooling->Redshifts[iz]);
-void construct_1d_table_from_3d(const struct part *restrict,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *, float *, int, float,
- int, int, float, int, int, double *, float *,
- float *);
+ get_redshift_index_previous = *z_index = iz;
-void construct_1d_print_table_from_3d_elements(
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *, float *, int,
- float, int, int, double *, float *, float *, float *);
+ break;
+ }
+ }
+ }
+}
-void construct_1d_table_from_3d_elements(
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *, float *, int,
- float, int, int, double *, float *, float *, float *);
+/*
+ * @brief Performs 1d interpolation (double precision)
+ *
+ * @param table Pointer to 1d table of values
+ * @param i Index of cell we are interpolating
+ * @param dx Offset within cell
+ */
+__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table,
+ int i, float dx) {
+
+ return (1 - dx) * table[i] + dx * table[i + 1];
+}
+
+/*
+ * @brief Performs 2d interpolation
+ *
+ * @param table Pointer to flattened 2d table of values
+ * @param i,j Indices of cell we are interpolating
+ * @param dx,dy Offset within cell
+ * @param nx,ny Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dy = 1 (used for calculating derivatives)
+ * @param lower Pointer to value set to the table value at
+ * the when dy = 0 (used for calculating derivatives)
+ */
+__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i,
+ int j, float dx,
+ float dy, int nx,
+ int ny, double *upper,
+ double *lower) {
+ float result;
+ int index[4];
-void construct_1d_table_from_4d(const struct part *restrict,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *, float *, int, float,
- int, int, float, int, int, float, int, int,
- double *, float *, float *);
+ index[0] = row_major_index_2d(i, j, nx, ny);
+ index[1] = row_major_index_2d(i, j + 1, nx, ny);
+ index[2] = row_major_index_2d(i + 1, j, nx, ny);
+ index[3] = row_major_index_2d(i + 1, j + 1, nx, ny);
-void construct_1d_print_table_from_4d_elements(
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *, float *, int,
- float, int, int, float, int, int, double *, float *, float *, float *);
+ 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]];
+ dy = 1.0;
+ *upper = (1 - dx) * (1 - dy) * table[index[0]] +
+ (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
+ dx * dy * table[index[3]];
+ dy = 0.0;
+ *lower = (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;
+}
+
+/*
+ * @brief Performs 3d interpolation
+ *
+ * @param table Pointer to flattened 3d table of values
+ * @param i,j,k Indices of cell we are interpolating
+ * @param dx,dy,dz Offset within cell
+ * @param nx,ny,nz Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dz = 1 (used for calculating derivatives)
+ * @param lower Pointer to value set to the table value at
+ * the when dz = 0 (used for calculating derivatives)
+ */
+__attribute__((always_inline)) INLINE float interpol_3d(
+ float *table, int i, int j, int k, float dx, float dy, float dz, int nx,
+ int ny, int nz, double *upper, double *lower) {
+ float result;
+ int index[8];
+
+ index[0] = row_major_index_3d(i, j, k, nx, ny, nz);
+ index[1] = row_major_index_3d(i, j, k + 1, nx, ny, nz);
+ index[2] = row_major_index_3d(i, j + 1, k, nx, ny, nz);
+ index[3] = row_major_index_3d(i, j + 1, k + 1, nx, ny, nz);
+ index[4] = row_major_index_3d(i + 1, j, k, nx, ny, nz);
+ index[5] = row_major_index_3d(i + 1, j, k + 1, nx, ny, nz);
+ index[6] = row_major_index_3d(i + 1, j + 1, k, nx, ny, nz);
+ index[7] = row_major_index_3d(i + 1, j + 1, k + 1, nx, ny, nz);
+
+ float table0 = table[index[0]];
+ float table1 = table[index[1]];
+ float table2 = table[index[2]];
+ float table3 = table[index[3]];
+ float table4 = table[index[4]];
+ float table5 = table[index[5]];
+ float table6 = table[index[6]];
+ float table7 = table[index[7]];
+
+ result = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+ (1 - dx) * (1 - dy) * dz * table1 +
+ (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+ dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+ dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+ dz = 1.0;
+ *upper = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+ (1 - dx) * (1 - dy) * dz * table1 +
+ (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+ dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+ dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+ dz = 0.0;
+ *lower = (1 - dx) * (1 - dy) * (1 - dz) * table0 +
+ (1 - dx) * (1 - dy) * dz * table1 +
+ (1 - dx) * dy * (1 - dz) * table2 + (1 - dx) * dy * dz * table3 +
+ dx * (1 - dy) * (1 - dz) * table4 + dx * (1 - dy) * dz * table5 +
+ dx * dy * (1 - dz) * table6 + dx * dy * dz * table7;
+
+ return result;
+}
+
+/*
+ * @brief Performs 4d interpolation
+ *
+ * @param table Pointer to flattened 4d table of values
+ * @param i,j,k,l Indices of cell we are interpolating
+ * @param dx,dy,dz,dw Offset within cell
+ * @param nx,ny,nz,nw Table dimensions
+ * @param upper Pointer to value set to the table value at
+ * the when dw = 1 when used for interpolating table
+ * depending on metal species, dz = 1 otherwise
+ * (used for calculating derivatives)
+ * @param upper Pointer to value set to the table value at
+ * the when dw = 0 when used for interpolating table
+ * depending on metal species, dz = 0 otherwise
+ * (used for calculating derivatives)
+ */
+__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 nx, int ny, int nz, int nw, double *upper, double *lower) {
+ float result;
+ int index[16];
+
+ index[0] = row_major_index_4d(i, j, k, l, nx, ny, nz, nw);
+ index[1] = row_major_index_4d(i, j, k, l + 1, nx, ny, nz, nw);
+ index[2] = row_major_index_4d(i, j, k + 1, l, nx, ny, nz, nw);
+ index[3] = row_major_index_4d(i, j, k + 1, l + 1, nx, ny, nz, nw);
+ index[4] = row_major_index_4d(i, j + 1, k, l, nx, ny, nz, nw);
+ index[5] = row_major_index_4d(i, j + 1, k, l + 1, nx, ny, nz, nw);
+ index[6] = row_major_index_4d(i, j + 1, k + 1, l, nx, ny, nz, nw);
+ index[7] = row_major_index_4d(i, j + 1, k + 1, l + 1, nx, ny, nz, nw);
+ index[8] = row_major_index_4d(i + 1, j, k, l, nx, ny, nz, nw);
+ index[9] = row_major_index_4d(i + 1, j, k, l + 1, nx, ny, nz, nw);
+ index[10] = row_major_index_4d(i + 1, j, k + 1, l, nx, ny, nz, nw);
+ index[11] = row_major_index_4d(i + 1, j, k + 1, l + 1, nx, ny, nz, nw);
+ index[12] = row_major_index_4d(i + 1, j + 1, k, l, nx, ny, nz, nw);
+ index[13] = row_major_index_4d(i + 1, j + 1, k, l + 1, nx, ny, nz, nw);
+ index[14] = row_major_index_4d(i + 1, j + 1, k + 1, l, nx, ny, nz, nw);
+ index[15] = row_major_index_4d(i + 1, j + 1, k + 1, l + 1, nx, ny, nz, nw);
+
+ float table0 = table[index[0]];
+ float table1 = table[index[1]];
+ float table2 = table[index[2]];
+ float table3 = table[index[3]];
+ float table4 = table[index[4]];
+ float table5 = table[index[5]];
+ float table6 = table[index[6]];
+ float table7 = table[index[7]];
+ float table8 = table[index[8]];
+ float table9 = table[index[9]];
+ float table10 = table[index[10]];
+ float table11 = table[index[11]];
+ float table12 = table[index[12]];
+ float table13 = table[index[13]];
+ float table14 = table[index[14]];
+ float table15 = table[index[15]];
+
+ result = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+ (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+ (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+ (1 - dx) * (1 - dy) * dz * dw * table3 +
+ (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+ (1 - dx) * dy * (1 - dz) * dw * table5 +
+ (1 - dx) * dy * dz * (1 - dw) * table6 +
+ (1 - dx) * dy * dz * dw * table7 +
+ dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+ dx * (1 - dy) * (1 - dz) * dw * table9 +
+ dx * (1 - dy) * dz * (1 - dw) * table10 +
+ dx * (1 - dy) * dz * dw * table11 +
+ dx * dy * (1 - dz) * (1 - dw) * table12 +
+ dx * dy * (1 - dz) * dw * table13 +
+ dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+ if (nw == 9) {
+ // interpolating metal species
+ dz = 1.0;
+ } else {
+ dw = 1.0;
+ }
+ *upper = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+ (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+ (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+ (1 - dx) * (1 - dy) * dz * dw * table3 +
+ (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+ (1 - dx) * dy * (1 - dz) * dw * table5 +
+ (1 - dx) * dy * dz * (1 - dw) * table6 +
+ (1 - dx) * dy * dz * dw * table7 +
+ dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+ dx * (1 - dy) * (1 - dz) * dw * table9 +
+ dx * (1 - dy) * dz * (1 - dw) * table10 +
+ dx * (1 - dy) * dz * dw * table11 +
+ dx * dy * (1 - dz) * (1 - dw) * table12 +
+ dx * dy * (1 - dz) * dw * table13 +
+ dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+ if (nw == 9) {
+ // interpolating metal species
+ dz = 0.0;
+ } else {
+ dw = 0.0;
+ }
+ *lower = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table0 +
+ (1 - dx) * (1 - dy) * (1 - dz) * dw * table1 +
+ (1 - dx) * (1 - dy) * dz * (1 - dw) * table2 +
+ (1 - dx) * (1 - dy) * dz * dw * table3 +
+ (1 - dx) * dy * (1 - dz) * (1 - dw) * table4 +
+ (1 - dx) * dy * (1 - dz) * dw * table5 +
+ (1 - dx) * dy * dz * (1 - dw) * table6 +
+ (1 - dx) * dy * dz * dw * table7 +
+ dx * (1 - dy) * (1 - dz) * (1 - dw) * table8 +
+ dx * (1 - dy) * (1 - dz) * dw * table9 +
+ dx * (1 - dy) * dz * (1 - dw) * table10 +
+ dx * (1 - dy) * dz * dw * table11 +
+ dx * dy * (1 - dz) * (1 - dw) * table12 +
+ dx * dy * (1 - dz) * dw * table13 +
+ dx * dy * dz * (1 - dw) * table14 + dx * dy * dz * dw * table15;
+
+ return result;
+}
+
+/*
+ * @brief Interpolates temperature from internal energy based on table and
+ * calculates the size of the internal energy cell for the specified
+ * internal energy.
+ *
+ * @param log_10_u Log base 10 of internal energy
+ * @param delta_u Pointer to size of internal energy cell
+ * @param z_i Redshift index
+ * @param n_h_i Hydrogen number density index
+ * @param He_i Helium fraction index
+ * @param d_z Redshift offset
+ * @param d_n_h Hydrogen number density offset
+ * @param d_He Helium fraction offset
+ * @param cooling Cooling data structure
+ * @param cosmo Cosmology data structure
+ */
+__attribute__((always_inline)) INLINE double eagle_convert_u_to_temp(
+ double log_10_u, float *delta_u, int z_i, int n_h_i, int He_i, float d_z,
+ float d_n_h, float d_He,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo) {
-void construct_1d_table_from_4d_elements(
- const struct part *restrict, const struct cooling_function_data *restrict,
- const struct cosmology *restrict, const struct phys_const *, float *, int,
- float, int, int, float, int, int, double *, float *, float *, float *);
+ int u_i;
+ float d_u, logT;
+ double upper, lower;
-double eagle_convert_temp_to_u_1d_table(
- double, float *, const struct cooling_function_data *restrict);
+ get_index_1d(cooling->Therm, cooling->N_Temp, log_10_u, &u_i, &d_u);
-double eagle_convert_u_to_temp(double, float *, int, int, int, float, float,
- float,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict);
+ if (cosmo->z > cooling->reionisation_redshift) {
+ logT = interpol_3d(cooling->table.photodissociation_cooling.temperature,
+ n_h_i, He_i, u_i, d_n_h, d_He, d_u, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp, &upper, &lower);
+ } else if (cosmo->z > cooling->Redshifts[cooling->N_Redshifts - 1]) {
+ logT = interpol_3d(cooling->table.no_compton_cooling.temperature, n_h_i,
+ He_i, u_i, d_n_h, d_He, d_u, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp, &upper, &lower);
+ } else {
+ logT = interpol_4d(cooling->table.element_cooling.temperature, z_i, n_h_i,
+ He_i, u_i, d_z, d_n_h, d_He, d_u, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp, &upper,
+ &lower);
+ }
-double eagle_convert_u_to_temp_1d_table(
- double, float *, double *, const struct cooling_function_data *restrict);
+ *delta_u =
+ pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
-void construct_1d_tables(int, float, int, float, int, float,
- const struct phys_const *restrict,
- const struct cosmology *restrict,
- const struct cooling_function_data *restrict,
- const struct part *restrict, float *, double *,
- double *, double *, double *, double *, double *,
- float *, float *);
+ return logT;
+}
-#endif /* SWIFT_INTERPOL_EAGLE_H */
+#endif
diff --git a/src/engine.c b/src/engine.c
index 51d43db964109ce2fc78f59a5a6e2ee2f1dad6b8..ce403b1e912553aa9b0238f1bf13c9e63b9033cb 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -4569,11 +4569,6 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
space_init_parts(s, e->verbose);
space_init_gparts(s, e->verbose);
- /* Update the cooling function */
- // if (e->policy & engine_policy_cooling)
- // cooling_update(e->physical_constants, e->internal_units, e->cosmology,
- // e->cooling_func);
-
/* Now, launch the calculation */
TIMER_TIC;
engine_launch(e);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 89e08edd130111d3e1e3fad4c1400b633a655ef6..d99b68f224f542dcdc60ae59fc6e2042ae20d9b7 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -20,7 +20,7 @@ AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS) $(GSL_INCS) $(FFTW_INCS)
AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
# List of programs and scripts to run in the test suite
-TESTS = testGreetings testCooling testMaths testReading.sh testSingle testKernel testSymmetry \
+TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
testActivePair.sh test27cells.sh test27cellsPerturbed.sh \
testParser.sh testSPHStep test125cells.sh test125cellsPerturbed.sh testFFT \
testAdiabaticIndex \
@@ -31,7 +31,7 @@ TESTS = testGreetings testCooling testMaths testReading.sh testSingle testKernel
testCbrt testCosmology testOutputList testFormat.sh
# List of test programs to compile
-check_PROGRAMS = testGreetings testCooling testReading testSingle testTimeIntegration \
+check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSPHStep testActivePair test27cells test27cells_subset test125cells testParser \
testKernel testFFT testInteractions testMaths \
testSymmetry testThreadpool \
@@ -47,8 +47,6 @@ $(check_PROGRAMS): ../src/.libs/libswiftsim.a
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
-testCooling_SOURCES = testCooling.c
-
testMaths_SOURCES = testMaths.c
testReading_SOURCES = testReading.c
diff --git a/tests/plots.py b/tests/plots.py
deleted file mode 100644
index 73b95bed5789a00409edb0c0b06f0577bbaf3086..0000000000000000000000000000000000000000
--- a/tests/plots.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import matplotlib.pyplot as plt
-
-elements = 11
-temperature = []
-cooling_rate = [[] for i in range(elements+1)]
-u = []
-fu = []
-length = 0
-file_in = open('cooling_output.dat', 'r')
-for line in file_in:
- data = line.split()
- temperature.append(float(data[0]))
- cooling_rate[0].append(-float(data[1]))
-
-file_in.close()
-
-for elem in range(elements):
- file_in = open('cooling_output_'+str(elem)+'.dat','r')
- for line in file_in:
- data = line.split()
- cooling_rate[elem+1].append(-float(data[0]))
- file_in.close()
-
-file_in = open('newton_output.dat', 'r')
-for line in file_in:
- data = line.split()
- u.append(float(data[0]))
- fu.append(float(data[1]))
-
-file_in.close()
-
-p0, = plt.plot(u,fu, color = 'k')
-#p1, = plt.plot(u,[0 for x in u], color = 'r')
-#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
-plt.xlim([1e13,2.0e14])
-plt.ylim([0e13,2e14])
-plt.xlabel('u')
-plt.ylabel('f(u)')
-plt.show()
-
-#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 = '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
deleted file mode 100644
index 40eb7cd4596da1fdec165b244745af3c6728e01a..0000000000000000000000000000000000000000
--- a/tests/testCooling.c
+++ /dev/null
@@ -1,232 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main() {
-
- /* Read the parameter file */
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_data chemistry_data;
- struct part p;
- struct xpart xp;
- struct phys_const internal_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
- enum table_index {
- EAGLE_Hydrogen = 0,
- EAGLE_Helium,
- EAGLE_Carbon,
- EAGLE_Nitrogen,
- EAGLE_Oxygen,
- EAGLE_Neon,
- EAGLE_Magnesium,
- EAGLE_Silicon,
- EAGLE_Iron
- };
-
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- /* And dump the parameters as used. */
- parser_write_params_to_file(params, "used_parameters.yml");
-
- units_init(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &internal_const);
-
- double hydrogen_number_density_cgs = 1e-1;
- double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
- temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
- u_cgs = 0;
- cooling_du_dt = 0;
- temperature_cgs = 0;
- newton_func = 0;
- // int n_t_i = 2000;
- dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
-
- char *output_filename = "cooling_output.dat";
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- char *output_filename2 = "temperature_output.dat";
- FILE *output_file2 = fopen(output_filename2, "w");
- if (output_file2 == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- char *output_filename3 = "newton_output.dat";
- FILE *output_file3 = fopen(output_filename3, "w");
- if (output_file3 == NULL) {
- printf("Error opening file!\n");
- 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);
-
- cosmology_init(params, &us, &internal_const, &cosmo);
- cosmology_print(&cosmo);
-
- 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));
-
- cooling_init(params, &us, &internal_const, &cooling);
- cooling_print(&cooling);
-
- // construct 1d table of cooling rates wrt temperature
- float H_plus_He_heat_table[176]; // WARNING sort out how it is
- // declared/allocated
- float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is
- // declared/allocated
- float temp_table[176]; // WARNING sort out how it is declared/allocated
- float element_cooling_table[9 * 176]; // WARNING sort out how it is
- // declared/allocated
- float element_electron_abundance_table[176]; // WARNING sort out how it is
- // declared/allocated
- construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_heating,
- H_plus_He_heat_table);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_electron_abundance,
- H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.temperature,
- temp_table);
- construct_1d_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, element_cooling_table);
- construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.electron_abundance,
- element_electron_abundance_table);
-
- //
- // printf("cooling table values \n");
- // for( int j=0; j < 176; j++) {
- // printf(" %.5e", H_plus_He_heat_table[j]+element_cooling_table[j]
-
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
- &internal_const) *
- cooling.number_density_scale;
- ratefact = inn_h * (XH / eagle_proton_mass_cgs);
- u_ini_cgs = pow(10.0, 14);
-
- // Compute contributions to cooling rate from different metals
- // for(int t_i = 0; t_i < n_t_i; t_i++){
- //
- // u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
- // u =
- // u_cgs/(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));
-
- // //cooling_du_dt =
- // eagle_print_metal_cooling_rate(&p,&cooling,&cosmo,&internal_const);
- // cooling_du_dt =
- // eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cooling,&cosmo,&internal_const);
- // temperature_cgs =
- // eagle_convert_u_to_temp(u_cgs,&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);
-
- // newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
- // fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
- // //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
- // \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
- //}
- // fclose(output_file);
- // fclose(output_file2);
- // fclose(output_file3);
-
- double dLambdaNet_du, LambdaNet; //, LambdaNext;
- float x_init, u_eq = 2.0e12;
- for (int j = 5; j < 6; j++) {
- float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
- temp_table, &p, &cooling,
- &cosmo, &internal_const),
- x, du;
- float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
- LambdaNet = eagle_cooling_rate_1d_table(
- u_ini, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
- &internal_const);
- float u_temp = u_ini + LambdaNet * ratefact * dt;
- /* RGB removed this **
- if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
- &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
- element_cooling_table, element_electron_abundance_table, temp_table, &p,
- &cooling, &cosmo, &internal_const);
- if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
- u_temp = u_ini;
- } else {
- u_temp =
- eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
- } */
- if (u_temp > u_eq) {
- x_init = log(u_temp);
- } else {
- x_init = log(u_eq);
- }
- x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, &p, &cosmo,
- &cooling, &internal_const, dt);
- printf(
- "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
- "final %.5e %.5e\n",
- u_ini, exp(x),
- eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
- &cosmo, &internal_const),
- eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
- &cosmo, &internal_const));
- }
-
- free(params);
-
- return 0;
-}
diff --git a/tests/testCoolingIntegration.c b/tests/testCoolingIntegration.c
deleted file mode 100644
index 925ca8a6f7889e8f29112514a538bf7fac8d0759..0000000000000000000000000000000000000000
--- a/tests/testCoolingIntegration.c
+++ /dev/null
@@ -1,265 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-
-#include "cooling.h"
-#include "hydro.h"
-#include "physical_constants.h"
-#include "swift.h"
-#include "units.h"
-
-int main() {
-
- /* Read the parameter file */
- struct swift_params *params = malloc(sizeof(struct swift_params));
- struct unit_system us;
- struct chemistry_data chemistry_data;
- struct part p;
- struct xpart xp;
- struct phys_const internal_const;
- struct cooling_function_data cooling;
- struct cosmology cosmo;
- // char *parametersFileName = "../examples/CoolingBox/coolingBox.yml";
- char *parametersFileName = "../examples/EAGLE_12/eagle_12.yml";
- enum table_index {
- EAGLE_Hydrogen = 0,
- EAGLE_Helium,
- EAGLE_Carbon,
- EAGLE_Nitrogen,
- EAGLE_Oxygen,
- EAGLE_Neon,
- EAGLE_Magnesium,
- EAGLE_Silicon,
- EAGLE_Iron
- };
-
- if (params == NULL) error("Error allocating memory for the parameter file.");
- message("Reading runtime parameters from file '%s'", parametersFileName);
- parser_read_file(parametersFileName, params);
-
- char output_filename[40];
- FILE **output_file = malloc(10 * sizeof(FILE *));
-
- /* And dump the parameters as used. */
- parser_write_params_to_file(params, "used_parameters.yml");
-
- units_init(&us, params, "InternalUnitSystem");
- phys_const_init(&us, params, &internal_const);
-
- double hydrogen_number_density_cgs = 1.582e-3;
- // double hydrogen_number_density_cgs = 1.0e-1;
- double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
- temperature_cgs, newton_func, ratefact;
- u_cgs = 0;
- cooling_du_dt = 0;
- temperature_cgs = 0;
- newton_func = 0;
- // int n_t_i = 2000;
-
- 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);
-
- cosmology_init(params, &us, &internal_const, &cosmo);
- cosmology_print(&cosmo);
-
- float u_ini = 3.357e15;
- u = u_ini / (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));
- cosmo.z = 0.0999744;
-
- cooling_init(params, &us, &internal_const, &cooling);
- cooling_print(&cooling);
-
- // construct 1d table of cooling rates wrt temperature
- float H_plus_He_heat_table[176]; // WARNING sort out how it is
- // declared/allocated
- float H_plus_He_electron_abundance_table[176]; // WARNING sort out how it is
- // declared/allocated
- float temp_table[176]; // WARNING sort out how it is declared/allocated
- float element_cooling_table[176]; // WARNING sort out how it is
- // declared/allocated
- float element_electron_abundance_table[176]; // WARNING sort out how it is
- // declared/allocated
- for (int k = 0; k < cooling.N_Temp; k++) H_plus_He_heat_table[k] = 0.0;
- for (int k = 0; k < cooling.N_Temp; k++)
- H_plus_He_electron_abundance_table[k] = 0.0;
- for (int k = 0; k < cooling.N_Temp; k++) temp_table[k] = 0.0;
- for (int k = 0; k < cooling.N_Temp; k++) element_cooling_table[k] = 0.0;
- for (int k = 0; k < cooling.N_Temp; k++)
- element_electron_abundance_table[k] = 0.0;
- // construct_1d_table_from_4d(&p,&cooling,&cosmo,&internal_const,cooling.table.element_cooling.H_plus_He_heating,H_plus_He_heat_table);
-
- float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
- // float inn_h =
- // chemistry_get_number_density(&p,&cosmo,chemistry_element_H,&internal_const)*cooling.number_density_scale;
- float inn_h = hydro_get_physical_density(&p, &cosmo) *
- units_cgs_conversion_factor(&us, UNIT_CONV_DENSITY) * XH /
- eagle_proton_mass_cgs;
- ratefact = inn_h * (XH / eagle_proton_mass_cgs);
- printf("XH inn_h ratefact %.5e %.5e %.5e\n", XH, inn_h, ratefact);
-
- double dLambdaNet_du, LambdaNet;
- float x_init;
- for (int j = 0; j < 1; j++) {
- // float u_ini =
- // eagle_convert_temp_to_u_1d_table(pow(10.0,0.5*(j+5)),temp_table,&p,&cooling,&cosmo,&internal_const);
- // float u_ini = 3.357e15;
- u = u_ini / (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));
-
- float x, du;
- double u_temp;
- // float dt = 2.0e-4*units_cgs_conversion_factor(&us,UNIT_CONV_TIME);
- float dt = 1.73798e-3 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
- construct_1d_table_from_4d_H_He(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_heating, H_plus_He_heat_table,
- &u_temp, u_ini, dt);
- construct_1d_table_from_4d(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.H_plus_He_electron_abundance,
- H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d(&p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.temperature,
- temp_table);
- construct_1d_table_from_4d_elements(
- &p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.metal_heating, element_cooling_table);
- construct_1d_table_from_3d(&p, &cooling, &cosmo, &internal_const,
- cooling.table.element_cooling.electron_abundance,
- element_electron_abundance_table);
-
- sprintf(output_filename, "%s%d%s", "cooling_integration_output_", j,
- ".dat");
- output_file[j] = fopen(output_filename, "w");
- if (output_file[j] == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
- for (int k = 0; k < cooling.N_Temp - 1; k++) {
- float lambda1 =
- H_plus_He_heat_table[k] + element_cooling_table[k] *
- H_plus_He_electron_abundance_table[k] /
- element_electron_abundance_table[k];
- float lambda2 = H_plus_He_heat_table[k + 1] +
- element_cooling_table[k + 1] *
- H_plus_He_electron_abundance_table[k + 1] /
- element_electron_abundance_table[k + 1];
- // printf("temperature %.5e, internal energy
- // %.5e\n",pow(10.0,temp_table[k]),
- // eagle_convert_temp_to_u_1d_table(pow(10.0,temp_table[k]),temp_table,&p,&cooling,&cosmo,&internal_const));
- float u2 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k + 1]),
- temp_table, &p, &cooling,
- &cosmo, &internal_const);
- float u1 = eagle_convert_temp_to_u_1d_table(pow(10.0, temp_table[k]),
- temp_table, &p, &cooling,
- &cosmo, &internal_const);
- float delta_u = u2 - u1;
- // fprintf(output_file[j], "%.5e %.5e %.5e\n", pow(10.0,temp_table[k]), u1
- // - u_ini - lambda1*ratefact*dt, (lambda2 - lambda1)/delta_u);
- fprintf(output_file[j], "%.5e %.5e %.5e\n", u1,
- u1 - u_ini - lambda1 * ratefact * dt,
- (lambda2 - lambda1) / delta_u);
- }
- fclose(output_file[j]);
-
- LambdaNet = eagle_cooling_rate_1d_table(
- u_ini, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
- &internal_const);
- double u_eq = 5.0e12;
- double u_temp_guess = u_ini + LambdaNet * ratefact * dt;
- printf(
- "u_guess, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du %.5e %.5e %.5e "
- "%.5e %.5e %.5e %.5e \n",
- u_temp, u_temp_guess, u_ini, LambdaNet, dLambdaNet_du, ratefact, dt);
-
- // if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
- // (LambdaNet >= 0 && u_temp > u_temp_guess))
- // u_temp = u_temp_guess;
- u_temp = u_temp_guess;
- if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
- u_temp = u_eq;
-
- x_init = log(u_temp);
- // int printflag = 1;
- // x =
- // newton_iter(x_init,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_table,element_electron_abundance_table,temp_table,&p,&cosmo,&cooling,&internal_const,dt,&printflag);
- x = bisection_iter(x_init, u_ini, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_cooling_table, element_electron_abundance_table,
- temp_table, &p, &cosmo, &cooling, &internal_const, dt);
- // printf("testing newton integration, u_ini, u %.5e %.5e, temperature
- // initial, final %.5e %.5e\n", u_ini, exp(x),
- // eagle_convert_u_to_temp_1d_table(u_ini,&du,temp_table,&p,&cooling,&cosmo,&internal_const),
- // eagle_convert_u_to_temp_1d_table(exp(x),&du,temp_table,&p,&cooling,&cosmo,&internal_const));
-
- u = u_ini;
- double u_next;
- int nt = 10000;
- float dt_sub = dt / nt;
- for (int t = 0; t < nt; t++) {
- LambdaNet = eagle_cooling_rate_1d_table(
- u, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
- &internal_const);
- u_next = u + LambdaNet * ratefact * dt_sub;
- printf(
- "here u_next u lambda_net ratefact dt_sub, du t %.5e %.5e %.5e %.5e "
- "%.5e %.5e %d\n",
- u_next, u, LambdaNet, ratefact, dt_sub, LambdaNet * ratefact * dt_sub,
- t);
- u = u_next;
- }
- printf(
- "testing newton integration, u_ini, u, u subcycle %.5e %.5e %.5e, "
- "temperature initial, final, subcycled %.5e %.5e %.5e\n",
- u_ini, exp(x), u,
- eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
- &cosmo, &internal_const),
- eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
- &cosmo, &internal_const),
- eagle_convert_u_to_temp_1d_table(u, &du, temp_table, &p, &cooling,
- &cosmo, &internal_const));
- }
-
- free(params);
-
- return 0;
-}