diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 0d83f718ced588a3ea6208ff33ee2937e172c28a..16c2fc55777eafb367896b4454f2c716538e6058 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/examples/CoolingBox/analytical_test.py
@@ -1,3 +1,4 @@
+import sys
import matplotlib
matplotlib.use("Agg")
from pylab import *
@@ -30,33 +31,23 @@ import numpy as np
import h5py as h5
import sys
-def interpol_lambda(temp_list,cooling_rate,temp):
- #print(temp,temp_list[0],temp_list[len(temp_list)-1])
- if temp < temp_list[0]:
+# function for interpolating 1d table of cooling rates
+def interpol_lambda(u_list,cooling_rate,u):
+ if u < u_list[0]:
return[cooling_rate[0]]
- if temp > temp_list[len(temp_list)-1]:
+ if u > u_list[len(u_list)-1]:
return[cooling_rate[len(cooling_rate)-1]]
j = 0
- while temp_list[j+1] < temp:
+ while u_list[j+1] < u:
j += 1
- cooling = cooling_rate[j]*(temp_list[j+1]-temp)/(temp_list[j+1]-temp_list[j]) + cooling_rate[j+1]*(temp-temp_list[j])/(temp_list[j+1]-temp_list[j])
+ cooling = cooling_rate[j]*(u_list[j+1]-u)/(u_list[j+1]-u_list[j]) + cooling_rate[j+1]*(u-u_list[j])/(u_list[j+1]-u_list[j])
return cooling
-def convert_u_to_temp_sol(u,rho):
- k_b = 1.38E-16 #boltzmann
- m_p = 1.67e-24 #proton mass
- gamma = 5.0/3.0
- n = 2.0*0.752*rho/m_p + 0.248*rho/(4.0*m_p)
- pressure = u*rho*(gamma - 1.0)
- temp = pressure/(k_b*n)
- return temp
-
# File containing the total energy
stats_filename = "./energy.txt"
# First snapshot
snap_filename = "coolingBox_0000.hdf5"
-snap_filename_mat = "../../../swiftsim_matthieu_cooling_v2/examples/CoolingBox/coolingBox_0000.hdf5"
# Some constants in cgs units
k_b = 1.38E-16 #boltzmann
@@ -67,9 +58,7 @@ T_init = 1.0e7
# Read the initial state of the gas
f = h5.File(snap_filename,'r')
-f_mat = h5.File(snap_filename_mat,'r')
rho = np.mean(f["/PartType0/Density"])
-rho_mat = np.mean(f["/PartType0/Density"])
pressure = np.mean(f["/PartType0/Pressure"])
# Read the units parameters from the snapshot
@@ -78,96 +67,86 @@ unit_mass = units.attrs["Unit mass in cgs (U_M)"]
unit_length = units.attrs["Unit length in cgs (U_L)"]
unit_time = units.attrs["Unit time in cgs (U_t)"]
unit_vel = unit_length/unit_time
-#hubble_param = 0.71
-hubble_param = 1.0
-
-unit_length = unit_length/hubble_param
-unit_mass = unit_mass/hubble_param
-
-yHe = 0.28
-temp_0 = 1.0e7
rho = rho*unit_mass/(unit_length**3)
-rho_mat = rho_mat*unit_mass/(unit_length**3)
# Read snapshots
if len(sys.argv) >= 4:
- nsnap = int(sys.argv[4])
+ nsnap = int(sys.argv[5])
else:
- nsnap = 501
+ print("Need to specify number of snapshots, defaulting to 100.")
+ nsnap = 100
npart = 4096
u_snapshots_cgs = zeros(nsnap)
u_part_snapshots_cgs = zeros((nsnap,npart))
t_snapshots_cgs = zeros(nsnap)
+scale_factor = zeros(nsnap)
for i in range(nsnap):
snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
- u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*unit_length**2/(unit_time**2)
+ u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*(unit_length**2)/(unit_time**2)
t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
-
-# Read snapshots
-#nsnap_mat = 25
-#npart = 4096
-#u_snapshots_cgs_mat = zeros(nsnap_mat)
-#u_part_snapshots_cgs_mat = zeros((nsnap_mat,npart))
-#t_snapshots_cgs_mat = zeros(nsnap_mat)
-#for i in range(nsnap_mat):
-# snap = h5.File("../../../swiftsim_matthieu_cooling_v2/examples/CoolingBox/coolingBox_%0.4d.hdf5"%i,'r')
-# u_part_snapshots_cgs_mat[i][:] = snap["/PartType0/InternalEnergy"][:]*unit_length**2/(unit_time**2)
-# t_snapshots_cgs_mat[i] = snap["/Header"].attrs["Time"] * unit_time
+ scale_factor[i] = snap["/Header"].attrs["Scale-factor"]
# calculate analytic solution. Since cooling rate is constant,
# temperature decrease in linear in time.
# read Lambda and temperatures from table
-temperature = []
+internal_energy = []
cooling_rate = []
file_in = open('cooling_output.dat', 'r')
for line in file_in:
data = line.split()
- temperature.append(float(data[0]))
+ internal_energy.append(float(data[0]))
cooling_rate.append(-float(data[1]))
tfinal = t_snapshots_cgs[nsnap-1]
-nt = 1e5
-dt = tfinal/nt
+tfirst = t_snapshots_cgs[0]
+nt = nsnap*10
+dt = (tfinal-tfirst)/nt
t_sol = np.zeros(int(nt))
-temp_sol = np.zeros(int(nt))
u_sol = np.zeros(int(nt))
lambda_sol = np.zeros(int(nt))
-u = np.mean(u_part_snapshots_cgs[0,:])
-temp_sol[0] = convert_u_to_temp_sol(u,rho)
-#print(u,temp_sol[0])
+u = np.mean(u_part_snapshots_cgs[0,:])/scale_factor[0]**2
u_sol[0] = u
+t_sol[0] = tfirst
+# integrate explicit ODE
for j in range(int(nt-1)):
t_sol[j+1] = t_sol[j] + dt
- Lambda_net = interpol_lambda(temperature,cooling_rate,u_sol[j])
- nH = 0.702*rho/(m_p)
- #nH = 1.0e-4
- u_next = u - Lambda_net*nH**2/rho*dt
- temp_sol[j+1] = convert_u_to_temp_sol(u_next,rho)
+ Lambda_net = interpol_lambda(internal_energy,cooling_rate,u_sol[j])
+ if j < 10:
+ print(u,Lambda_net)
+ if int(sys.argv[4]) == 1:
+ nH = 0.702*rho/(m_p)/scale_factor[0]**3
+ ratefact = nH*0.702/m_p
+ else:
+ nH = 0.752*rho/(m_p)/scale_factor[0]**3
+ ratefact = nH*0.752/m_p
+ u_next = u - Lambda_net*ratefact*dt
u_sol[j+1] = u_next
- #lambda_sol[j] = Lambda_net
u = u_next
-
-mean_temp = np.zeros(nsnap)
+u_sol = u_sol*scale_factor[0]**2
+
+# swift solution
mean_u = np.zeros(nsnap)
for j in range(nsnap):
- mean_temp[j] = convert_u_to_temp_sol(np.mean(u_part_snapshots_cgs[j,:]),rho)
mean_u[j] = np.mean(u_part_snapshots_cgs[j,:])
-#mean_temp_mat = np.zeros(nsnap_mat)
-#for j in range(nsnap_mat):
-# mean_temp_mat[j] = convert_u_to_temp_sol(np.mean(u_part_snapshots_cgs_mat[j,:]),rho_mat)
+# plot and save results
+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')
+p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean alexei')
l = legend(handles = [p1,p2])
xlabel("${\\rm{Time~[s]}}$", labelpad=0, fontsize = 14)
ylabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
-title('$n_h = 1 \\rm{cm}^{-3}, z = 0$, zero metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
-#plt.ylim([8.0e11,3.0e12])
+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)
+ 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)
+ name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_zero_metal.png"
-savefig("energy.png", dpi=200)
+savefig(name, dpi=200)
print('Final internal energy ode, swift, error ' + str(u_sol[int(nt)-1]) + ' ' + str(mean_u[nsnap-1]) + ' ' + str( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])))
-#print(u_sol[int(nt)-1], mean_u[nsnap-1], (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1]))
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index aa0597c05fecabf007af0b4cda84bbf4f324c2b0..03f5626023bf1570caa65df52673c8b5565ec4a3 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 = 3.271e7 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 4.5e9 # Pressure in code units (at 10^6K)
+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)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 0d3959790d77dcf4512035eb95348a45fe68121b..4c97e82d3e9019869f6517d535dcb844f9660f43 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -21,7 +21,7 @@ then
fi
# Run SWIFT
-../swift -s -c -C -t 16 -n 2000 coolingBox.yml
+../swift -s -c -C -t 16 -n 1000 coolingBox.yml
# Check energy conservation and cooling rate
# python energy_plot.py
diff --git a/examples/CoolingBox/test_script.sh b/examples/CoolingBox/test_script.sh
index e3673256289bf7dd7fe16d493b5f4e478122e724..31bddb8c406397081e7913785785a482f82b09c2 100644
--- a/examples/CoolingBox/test_script.sh
+++ b/examples/CoolingBox/test_script.sh
@@ -5,13 +5,13 @@ max_number() {
}
# loop over redshift
-for z in 0.1; do
+for z in $(seq 0.01 5.0 25.01); do
# loop over solar and zero metal abundances
- for solar in 1; do
+ for solar in {0..1}; do
# loop over log_10 hydrogen number density
- for nh in -1; do
+ for nh_exp in $(seq -3 2.0 -1); do
#change parameters in yml file for calculating explicit ode solution of cooling
- cd ../CoolingTest_Alexei
+ cd ../CoolingRates
if [ $solar == 1 ]
then
sed -i "/InitMetallicity: / s/: \+[[:alnum:],\.,-]\+/: 0.014/g" testCooling.yml
@@ -41,10 +41,10 @@ for z in 0.1; do
sed -i "/CalciumOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.0941736/g" testCooling.yml
sed -i "/SulphurOverSilicon: / s/: \+[[:alnum:],\.,-]\+/: 0.6054160/g" testCooling.yml
fi
- # ./testCooling metals $nh
- ./testCooling -z $z -d $nh -t
+ rm cooling_*.dat
+ ./testCooling -z $z -d $nh_exp
cd ../CoolingBox
- cp ../CoolingTest_Alexei/cooling_output.dat ./
+ cp ../CoolingRates/cooling_output.dat ./
#change parameters in coolingBox.yml
if [ $solar == 1 ]
@@ -79,7 +79,7 @@ for z in 0.1; do
# set starting, ending redshift, how often to write to file
a_begin=$(python -c "print 1.0/(1.0+$z)")
- a_end=$(python -c "print min(1.0, $a_begin*1.001)")
+ a_end=$(python -c "print min(1.0, $a_begin*1.0001)")
first_ouput_a=$a_begin
delta_a=$(python -c "print 1.0 + ($a_end/$a_begin - 1.0)/500.")
sed -i "/a_begin: / s/: \+[[:alnum:],\.,-]\+/: $a_begin/g" coolingBox.yml
@@ -88,26 +88,28 @@ for z in 0.1; do
sed -i "/delta_time: / s/: \+[[:alnum:],\.,-]\+/: $delta_a/g" coolingBox.yml
# change hydrogen number density
- sed -i "/^rho =/ s/= \S*/= 3.2e$(expr $nh + 7)/g" makeIC.py
- for pressure in 7; do
+ nh=$(python -c "print 3.555*10.0**($nh_exp+7)/(1.0+$z)**3")
+ sed -i "/^rho =/ s/= \S*/= $nh/g" makeIC.py
+ for pressure_index in 7; do
+ pressure=$(python -c "print 4.5*10.0**($pressure_index + $nh_exp + 3)/(1.0+$z)")
# change pressure (and hence energy)
- sed -i "/^P =/ s/= \S*/= 4.5e$(expr $pressure + $nh + 2)/g" makeIC.py
+ sed -i "/^P =/ s/= \S*/= $pressure/g" makeIC.py
python makeIC.py
+ rm coolingBox_*hdf5
# run cooling box
./run.sh
max=0
- for file in $( ls coolingBox_*.hdf5 ); do
+ for file in $( ls -lth coolingBox_*.hdf5 | head -n 1 ); do
f=${file//hdf5/}
f=${f//[!0-9]/}
- max=$(max_number $max $f)
done
- frames=$((10#$max))
+ frames=$((10#$f))
echo "number of frames $frames"
# check if everything worked and create plots
- python analytical_test.py $nh $pressure $solar $frames
+ python analytical_test.py $z $nh_exp $pressure_index $solar $frames
done
done
done
diff --git a/examples/CoolingRates/plots.py b/examples/CoolingRates/plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..969e6de713a4c421cc14d5f88ef68bafaab59a41
--- /dev/null
+++ b/examples/CoolingRates/plots.py
@@ -0,0 +1,58 @@
+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/CoolingRates/plots_nh.py b/examples/CoolingRates/plots_nh.py
new file mode 100644
index 0000000000000000000000000000000000000000..ab3300750e2497ec8026dbe2d5c6c7d358f87831
--- /dev/null
+++ b/examples/CoolingRates/plots_nh.py
@@ -0,0 +1,47 @@
+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/CoolingRates/testCooling.c b/examples/CoolingRates/testCooling.c
index 22ca23423ebc87cf8ca5a9ba0b08762641df13ab..4f141097a99d1b73e4eac4c6ffb3f3098cebdd9f 100644
--- a/examples/CoolingRates/testCooling.c
+++ b/examples/CoolingRates/testCooling.c
@@ -69,6 +69,299 @@ 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,
@@ -154,6 +447,14 @@ 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] /
@@ -163,48 +464,139 @@ 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);
- int nt = 250;
+ 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
- // 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);
- }
+ //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));
- 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");
+ // 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);
+ //}
+
+ // 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/CoolingRates/testCooling.yml b/examples/CoolingRates/testCooling.yml
index 3010636433d35421009539d1858722ea89820bc9..db3a7cc00a529dcc403d4c24c11b30186baa7896 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.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
+ 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
CalciumOverSilicon: 0.0941736
SulphurOverSilicon: 0.6054160
diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h
index 5d7e0e8883d81996c1f8e10288f7d63bdeb9c365..04bc6371973743bfdd88c68bde38a2d95ec0fe80 100644
--- a/src/chemistry/EAGLE/chemistry.h
+++ b/src/chemistry/EAGLE/chemistry.h
@@ -214,8 +214,6 @@ chemistry_get_number_density(const struct part* restrict p,
double element_mass = internal_const->const_proton_mass * atomic_number;
number_density = p->chemistry_data.metal_mass_fraction[elem] *
hydro_get_physical_density(p, cosmo) / element_mass;
- // number_density =
- // p->chemistry_data.metal_mass_fraction[elem]*hydro_get_comoving_density(p)/element_mass;
return number_density;
}
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 419a2fb162926afe56fe6e0473a80300bc3c7050..0f0b90b17140ed38570e09733318c53f6a7c6abc 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -17,8 +17,8 @@
*
******************************************************************************/
/**
- * @file src/cooling/EAGLE/cooling.h
- * @brief EAGLE cooling function
+ * @file src/cooling/EAGLE/cooling.c
+ * @brief EAGLE cooling functions
*/
/* Config parameters. */
@@ -45,6 +45,9 @@
#include "units.h"
const float explicit_tolerance = 0.05;
+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
/*
* @brief calculates heating due to helium reionization
@@ -116,7 +119,7 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
double T_gam, solar_electron_abundance, solar_electron_abundance1, solar_electron_abundance2, elem_cool1, elem_cool2;
double n_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
internal_const) *
- cooling->number_density_scale; // chemistry_data
+ cooling->number_density_scale;
double z = cosmo->z;
double cooling_rate = 0.0, temp_lambda, temp_lambda1, temp_lambda2;
float du;
@@ -135,8 +138,8 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
// the temperature table.
temp = eagle_convert_u_to_temp(log_10_u, &du, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, p,
- cooling, cosmo, internal_const);
+ dz, d_n_h, d_He,
+ cooling, cosmo);
get_index_1d(cooling->Temp, cooling->N_Temp, temp, &temp_i, &d_temp);
/* ------------------ */
@@ -305,14 +308,12 @@ __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, int z_index,
- float dz, int n_h_i, float d_n_h, int He_i, float d_He,
+ 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,
- float *solar_ratio) {
+ double *element_lambda) {
double T_gam, solar_electron_abundance;
double n_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
internal_const) *
@@ -330,8 +331,7 @@ eagle_metal_cooling_rate_1d_table(
// 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, p, cooling, cosmo,
- internal_const);
+ 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);
/* ------------------ */
@@ -480,15 +480,14 @@ __attribute__((always_inline)) INLINE double eagle_cooling_rate_1d_table(
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, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
+ 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.
- * Prints contribution from each element to cooling rate for testing purposes.
+ * Writes to file 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.
*
@@ -554,12 +553,6 @@ eagle_print_metal_cooling_rate(
* @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
@@ -570,12 +563,10 @@ __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, int z_index, float dz, int n_h_i, float d_n_h, int He_i,
- float d_He, const struct part *restrict p,
+ 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,
- float *abundance_ratio) {
+ 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) *
@@ -597,8 +588,7 @@ eagle_print_metal_cooling_rate_1d_table(
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, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
+ 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]);
}
@@ -614,13 +604,6 @@ eagle_print_metal_cooling_rate_1d_table(
*
* @param logu_init Initial guess for log(internal energy)
* @param u_ini Internal energy at beginning of hydro step
- * @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
@@ -634,18 +617,15 @@ eagle_print_metal_cooling_rate_1d_table(
* @param dt Hydro timestep
*/
__attribute__((always_inline)) INLINE float bisection_iter(
- float logu_init, double u_ini, 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 logu_init, double u_ini, int z_index,
float dz, int n_h_i, float d_n_h, int He_i, float d_He, float He_reion_heat,
struct part *restrict p, const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct phys_const *restrict phys_const,
- float *abundance_ratio, float dt, float *ub, float *lb) {
- double logu_upper, logu_lower, logu_next, LambdaNet_upper, LambdaNet_lower, LambdaNet_next, LambdaNet, dLambdaNet_du;
+ float *abundance_ratio, float dt) {
+ double u_upper, u_lower, u_next, LambdaNet, dLambdaNet_du;
int i = 0;
- const float bisection_tolerance = 1.0e-8;
- const double bracket_factor = log(sqrt(1.1));
+ double u_init = exp(logu_init);
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
@@ -659,64 +639,60 @@ __attribute__((always_inline)) INLINE float bisection_iter(
double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
// Bracketing
- logu_lower = logu_init;
- logu_upper = logu_init;
+ u_lower = u_init;
+ u_upper = u_init;
LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_init, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ log(u_init), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
if (LambdaNet < 0){
- logu_lower -= bracket_factor;
- logu_upper += bracket_factor;
+ u_lower /= bracket_factor;
+ u_upper *= bracket_factor;
- while(LambdaNet < 0) {
- logu_lower -= bracket_factor;
- logu_upper -= bracket_factor;
+ LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
+ log(u_lower), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
+ while(u_lower - u_ini - LambdaNet*ratefact*dt > 0) {
+ u_lower /= bracket_factor;
+ u_upper /= bracket_factor;
LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_lower, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ log(u_lower), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
}
- LambdaNet_lower = LambdaNet;
- LambdaNet_upper = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_upper, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
} else {
- logu_lower -= bracket_factor;
- logu_upper += bracket_factor;
+ u_lower /= bracket_factor;
+ u_upper *= bracket_factor;
- while(LambdaNet > 0) {
- logu_lower += bracket_factor;
- logu_upper += bracket_factor;
+ LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
+ log(u_upper), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
+ while(u_upper - u_ini - LambdaNet*ratefact*dt < 0) {
+ u_lower *= bracket_factor;
+ u_upper *= bracket_factor;
LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_upper, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ log(u_upper), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
}
- LambdaNet_upper = LambdaNet;
- LambdaNet_lower = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_lower, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
}
// bisection iterations
i = 0;
do {
- logu_next = 0.5 * (logu_lower + logu_upper);
- LambdaNet_next = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
- logu_next, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ u_next = 0.5 * (u_lower + u_upper);
+ LambdaNet = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate(
+ log(u_next), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- if (LambdaNet_lower * LambdaNet_next < 0) {
- logu_upper = logu_next;
- LambdaNet_upper = LambdaNet_next;
+ if (u_next - u_ini - LambdaNet*ratefact*dt > 0.0) {
+ u_upper = u_next;
} else {
- logu_lower = logu_next;
- LambdaNet_lower = LambdaNet_next;
+ u_lower = u_next;
}
i++;
- } while (fabs(logu_lower - logu_upper) > bisection_tolerance && i < 50*eagle_max_iterations);
+ } while (fabs(u_upper - u_lower)/u_next > bisection_tolerance && i < 50*eagle_max_iterations);
- if (i >= 50*eagle_max_iterations) n_eagle_cooling_rate_calls_4++;
+ if (i >= 50*eagle_max_iterations) error("Particle id %llu failed to converge", p->id);
- return logu_upper;
+ return log(u_upper);
}
@@ -740,9 +716,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
* @param phys_const Physical constants structure
* @param abundance_ratio Array of ratios of metal abundance to solar
* @param dt timestep
- * @param printflag Flag to identify if scheme failed to converge
- * @param lb lower bound to be used by bisection scheme
- * @param ub upper bound to be used by bisection scheme
+ * @param bisection_flag Flag to identify if scheme failed to converge
*/
__attribute__((always_inline)) INLINE float newton_iter(
float logu_init, double u_ini, int z_index,
@@ -751,12 +725,10 @@ __attribute__((always_inline)) INLINE float newton_iter(
const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct phys_const *restrict phys_const,
- float *abundance_ratio, float dt, int *printflag,
- float *lb, float *ub) {
+ float *abundance_ratio, float dt, int *bisection_flag) {
/* this routine does the iteration scheme, call one and it iterates to
* convergence */
- const float newton_tolerance = 1.0e-2;
double logu, logu_old;
double dLambdaNet_du, LambdaNet;
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
@@ -773,30 +745,17 @@ __attribute__((always_inline)) INLINE float newton_iter(
logu_old = logu_init;
logu = logu_old;
int i = 0;
- // float log_table_bound_high = 43.749, log_table_bound_low = 20.723;
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;
do /* iterate to convergence */
{
- // Count how many steps
- n_eagle_cooling_rate_calls_1++;
-
logu_old = logu;
LambdaNet = (He_reion_heat / (dt * ratefact)) +
eagle_cooling_rate(
logu_old, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- // Assign to upper, lower bounds the highest or lowest values of
- // cooling rate seen so far.
- if (LambdaNet < 0) {
- *ub = fmin(*ub,logu_old);
- } else if (LambdaNet > 0) {
- *lb = fmax(*lb,logu_old);
- }
-
-
// Newton iteration.
logu = logu_old - (1.0 - u_ini * exp(-logu_old) -
LambdaNet * ratefact * dt * exp(-logu_old)) /
@@ -814,10 +773,8 @@ __attribute__((always_inline)) INLINE float newton_iter(
i++;
} while (fabs(logu - logu_old) > newton_tolerance && i < eagle_max_iterations);
if (i >= eagle_max_iterations) {
- // count how many failed to converge
- n_eagle_cooling_rate_calls_3++;
// flag to trigger bisection scheme
- if (*printflag == 0) *printflag = 1;
+ *bisection_flag = 1;
}
return logu;
@@ -883,10 +840,10 @@ void cooling_cool_part(
double ratefact, u, LambdaNet, LambdaTune = 0, dLambdaNet_du;
int He_i, n_h_i;
float d_He, d_n_h;
- const float hydro_du_dt = hydro_get_internal_energy_dt(p);
+ const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
double u_old = (hydro_get_physical_internal_energy(p, cosmo)
- + hydro_get_internal_energy_dt(p) * dt) *
+ + hydro_du_dt * dt) *
cooling->internal_energy_scale;
get_redshift_index(cosmo->z, &z_index, &dz, cooling);
@@ -909,6 +866,7 @@ void cooling_cool_part(
/* 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);
+
/* set helium and hydrogen reheating term */
if (cooling->he_reion == 1) LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
@@ -922,51 +880,24 @@ void cooling_cool_part(
log(u_ini), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- if (p->id == 1823) printf("Eagle cooling.h id dt dt normalised %llu %.5e %.5e\n", p->id, dt, dt/units_cgs_conversion_factor(us, UNIT_CONV_TIME));
-
if (fabs(ratefact * LambdaNet * dt) < explicit_tolerance * u_old) {
// if cooling rate is small, take the explicit solution
u = u_old + ratefact * LambdaNet * dt;
} else {
- // count how many times we use newton scheme
- n_eagle_cooling_rate_calls_2++;
- float logu, lower_bound = cooling->Therm[0]/eagle_log_10_e, upper_bound = cooling->Therm[cooling->N_Temp-1]/eagle_log_10_e;
+ float logu;
double u_eq = u_ini;
if (dt > 0) {
+ // newton method
double log_u_temp = log(u_eq);
-
- int printflag = 0;
+ int bisection_flag = 0;
logu = newton_iter(log_u_temp, u_ini, z_index, dz,
n_h_i, d_n_h, He_i, d_He, LambdaTune, p, cosmo, cooling, phys_const,
- abundance_ratio, dt, &printflag, &lower_bound, &upper_bound);
- if (printflag == 1) {
+ abundance_ratio, dt, &bisection_flag);
+ if (bisection_flag == 1) {
// newton method didn't work, so revert to bisection
-
- // construct 1d table of cooling rates wrt temperature
- double temp_table[cooling->N_Temp]; // WARNING sort out how it is declared/allocated
- double H_plus_He_heat_table[cooling->N_Temp];
- double H_plus_He_electron_abundance_table[cooling->N_Temp];
- double element_cooling_table[cooling->N_Temp];
- double element_cooling_print_table[cooling->N_Elements*cooling->N_Temp];
- double element_electron_abundance_table[cooling->N_Temp];
- construct_1d_tables(z_index, dz, n_h_i, d_n_h, He_i, d_He,
- phys_const, us, cosmo, cooling, p,
- abundance_ratio,
- H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_cooling_table,
- element_cooling_print_table,
- element_electron_abundance_table,
- temp_table,
- &upper_bound,
- &lower_bound);
- // perform bisection scheme
- logu = bisection_iter(log_u_temp,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_print_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,LambdaTune,p,cosmo,cooling,phys_const,abundance_ratio,dt,&upper_bound,&lower_bound);
-
- // do not trigger bisection again.
- printflag = 2;
+ logu = bisection_iter(log_u_temp,u_ini,z_index,dz,n_h_i,d_n_h,He_i,d_He,LambdaTune,p,cosmo,cooling,phys_const,abundance_ratio,dt);
}
u = exp(logu);
}
@@ -975,11 +906,12 @@ void cooling_cool_part(
// calculate du/dt
float cooling_du_dt = 0.0;
if (dt > 0) {
- cooling_du_dt = (u - u_ini) / dt / cooling->power_scale;
+ cooling_du_dt = (u - u_ini) / dt / cooling->power_scale / cosmo->a2_inv;
}
/* Update the internal energy time derivative */
- hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
+
+ hydro_set_physical_internal_energy_dt(p, cosmo, hydro_du_dt + cooling_du_dt);
/* Store the radiated energy */
xp->cooling_data.radiated_energy +=
@@ -1123,7 +1055,7 @@ INLINE void cooling_update(const struct phys_const* phys_const,
* @param cooling Cooling data structure containing properties to initialize
*/
INLINE void cooling_init_backend(
- const struct swift_params *parameter_file, const struct unit_system *us,
+ struct swift_params *parameter_file, const struct unit_system *us,
const struct phys_const *phys_const,
struct cooling_function_data *cooling) {
@@ -1152,7 +1084,6 @@ INLINE void cooling_init_backend(
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname, cooling);
cooling->table = eagle_readtable(cooling->cooling_table_path, cooling);
- printf("Eagle cooling.h read table \n");
// allocate array for solar abundances
cooling->solar_abundances = malloc(cooling->N_Elements * sizeof(double));
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index 1266b1aede587274bbdd8750e5c09b09bcbf8bb0..c6d4cb269ca50e34834453fdf04a6464fc86eb9d 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -21,39 +21,17 @@
/**
* @file src/cooling/EAGLE/cooling.h
- * @brief EAGLE cooling function
+ * @brief EAGLE cooling function declarations
*/
-/* Config parameters. */
-#include "../config.h"
-
-/* Some standard headers. */
-#include <float.h>
-#include <hdf5.h>
-#include <math.h>
-#include <time.h>
-
/* Local includes. */
-#include "chemistry.h"
#include "cooling_struct.h"
+#include "cosmology.h"
#include "eagle_cool_tables.h"
-#include "error.h"
-#include "hydro.h"
-#include "io_properties.h"
-#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"
-/* number of calls to eagle cooling rate */
-extern int n_eagle_cooling_rate_calls_1;
-extern int n_eagle_cooling_rate_calls_2;
-extern int n_eagle_cooling_rate_calls_3;
-extern int n_eagle_cooling_rate_calls_4;
-
-static int get_redshift_index_first_call = 0;
-static int get_redshift_index_previous = -1;
-
enum hdf5_allowed_types {
hdf5_short,
hdf5_int,
@@ -63,2248 +41,126 @@ enum hdf5_allowed_types {
hdf5_char
};
+double eagle_helium_reionization_extraheat(double, double, const struct cooling_function_data *restrict);
-/**
- * @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) {
- int index = i * ny + j;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
-#endif
- return index;
-}
-
-/**
- * @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) {
- int index = i * ny * nz + j * nz + k;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
- assert(k < nz);
-#endif
- return index;
-}
-
-/**
- * @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) {
- int index = i * ny * nz * nw + j * nz * nw + k * nw + l;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
- assert(k < nz);
- assert(l < nw);
-#endif
- return index;
-}
-
-/*
- * @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) {
- float dxm1;
- const float EPS = 1.e-4;
-
- dxm1 = (float)(ntable - 1) / (table[ntable - 1] - table[0]);
-
- if ((float)x <= table[0] + EPS) {
- *i = 0;
- *dx = 0;
- } else if ((float)x >= table[ntable - 1] - EPS) {
- *i = ntable - 2;
- *dx = 1;
- } else {
- *i = (int)floor(((float)x - table[0]) * dxm1);
- *dx = ((float)x - table[*i]) * dxm1;
- }
-}
-
-/*
- * @brief This routine returns the position i of a value z in the redshift table
- * and the
- * displacement dx needed for the interpolation.
- *
- * @param z Redshift whose position within the redshift array we are interested
- * in
- * @param z_index i Pointer to the index whose corresponding redshift
- * is the greatest value in the redshift table less than x
- * @param dx Pointer to offset of z within redshift cell
- * @param cooling Pointer to cooling structure (contains redshift table)
- */
-__attribute__((always_inline)) INLINE static void get_redshift_index(
- float z, int *z_index, float *dz,
- const struct cooling_function_data *restrict cooling) {
- int i, iz;
-
- if (get_redshift_index_first_call == 0) {
- get_redshift_index_first_call = 1;
- get_redshift_index_previous = cooling->N_Redshifts - 2;
-
- /* 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");
- }
- }
-
- /* 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]);
-
- get_redshift_index_previous = *z_index = iz;
-
- break;
- }
- }
- }
-}
-
-/*
- * @brief Performs 1d interpolation
- *
- * @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 float interpol_1d(float *table, int i,
- float dx) {
- float result;
-
- result = (1 - dx) * table[i] + dx * table[i + 1];
-
- return result;
-}
-
-/*
- * @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) {
- double result;
-
- result = (1 - dx) * table[i] + dx * table[i + 1];
-
- return result;
-}
-
-/*
- * @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
- */
-__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i,
- int j, float dx,
- float dy, int nx,
- int ny) {
- float result;
- int index[4];
-
- 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);
-#ifdef SWIFT_DEBUG_CHECKS
- if (index[0] >= nx * ny || index[0] < 0)
- fprintf(stderr, "index 0 out of bounds %d, i,j %d, %d, table size %d\n",
- index[0], i, j, nx * ny);
- if (index[1] >= nx * ny || index[1] < 0)
- fprintf(stderr, "index 1 out of bounds %d, i,j %d, %d, table size %d\n",
- index[1], i, j + 1, nx * ny);
- if (index[2] >= nx * ny || index[2] < 0)
- fprintf(stderr, "index 2 out of bounds %d, i,j %d, %d, table size %d\n",
- index[2], i + 1, j, nx * ny);
- if (index[3] >= nx * ny || index[3] < 0)
- fprintf(stderr, "index 3 out of bounds %d, i,j %d, %d, table size %d\n",
- index[3], i + 1, j + 1, nx * ny);
-#endif
-
- result = (1 - dx) * (1 - dy) * table[index[0]] +
- (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
- dx * dy * table[index[3]];
-
- return result;
-}
-
-/*
- * @brief Performs 2d interpolation (double precision)
- *
- * @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
- */
-__attribute__((always_inline)) INLINE double interpol_2d_dbl(
- double *table, int i, int j, double dx, double dy, int nx, int ny) {
- double result;
- int index[4];
-
- 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);
-#ifdef SWIFT_DEBUG_CHECKS
- if (index[0] >= nx * ny || index[0] < 0)
- fprintf(stderr, "index 0 out of bounds %d, i,j %d, %d, table size %d\n",
- index[0], i, j, nx * ny);
- if (index[1] >= nx * ny || index[1] < 0)
- fprintf(stderr, "index 1 out of bounds %d, i,j %d, %d, table size %d\n",
- index[1], i, j + 1, nx * ny);
- if (index[2] >= nx * ny || index[2] < 0)
- fprintf(stderr, "index 2 out of bounds %d, i,j %d, %d, table size %d\n",
- index[2], i + 1, j, nx * ny);
- if (index[3] >= nx * ny || index[3] < 0)
- fprintf(stderr, "index 3 out of bounds %d, i,j %d, %d, table size %d\n",
- index[3], i + 1, j + 1, nx * ny);
-#endif
-
- result = (1 - dx) * (1 - dy) * table[index[0]] +
- (1 - dx) * dy * table[index[1]] + dx * (1 - dy) * table[index[2]] +
- dx * dy * table[index[3]];
-
- return result;
-}
-
-/*
- * @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
- */
-__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);
-#ifdef SWIFT_DEBUG_CHECKS
- if (index[0] >= nx * ny * nz || index[0] < 0)
- fprintf(stderr,
- "index 0 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[0], i, j, k, nx * ny * nz);
- if (index[1] >= nx * ny * nz || index[1] < 0)
- fprintf(stderr,
- "index 1 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[1], i, j, k + 1, nx * ny * nz);
- if (index[2] >= nx * ny * nz || index[2] < 0)
- fprintf(stderr,
- "index 2 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[2], i, j + 1, k, nx * ny * nz);
- if (index[3] >= nx * ny * nz || index[3] < 0)
- fprintf(stderr,
- "index 3 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[3], i, j + 1, k + 1, nx * ny * nz);
- if (index[4] >= nx * ny * nz || index[4] < 0)
- fprintf(stderr,
- "index 4 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[4], i + 1, j, k, nx * ny * nz);
- if (index[5] >= nx * ny * nz || index[5] < 0)
- fprintf(stderr,
- "index 5 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[5], i + 1, j, k + 1, nx * ny * nz);
- if (index[6] >= nx * ny * nz || index[6] < 0)
- fprintf(stderr,
- "index 6 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[6], i + 1, j + 1, k, nx * ny * nz);
- if (index[7] >= nx * ny * nz || index[7] < 0)
- fprintf(stderr,
- "index 7 out of bounds %d, i,j,k %d, %d, %d, table size %d\n",
- index[7], i + 1, j + 1, k + 1, nx * ny * nz);
-#endif
-
- 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
- */
-__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);
-#ifdef SWIFT_DEBUG_CHECKS
- if (index[0] >= nx * ny * nz * nw || index[0] < 0)
- fprintf(
- stderr,
- "index 0 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[0], i, j, k, l, nx * ny * nz * nw);
- if (index[1] >= nx * ny * nz * nw || index[1] < 0)
- fprintf(
- stderr,
- "index 1 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[1], i, j, k, l + 1, nx * ny * nz * nw);
- if (index[2] >= nx * ny * nz * nw || index[2] < 0)
- fprintf(
- stderr,
- "index 2 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[2], i, j, k + 1, l, nx * ny * nz * nw);
- if (index[3] >= nx * ny * nz * nw || index[3] < 0)
- fprintf(
- stderr,
- "index 3 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[3], i, j, k + 1, l + 1, nx * ny * nz * nw);
- if (index[4] >= nx * ny * nz * nw || index[4] < 0)
- fprintf(
- stderr,
- "index 4 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[4], i, j + 1, k, l, nx * ny * nz * nw);
- if (index[5] >= nx * ny * nz * nw || index[5] < 0)
- fprintf(
- stderr,
- "index 5 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[5], i, j + 1, k, l + 1, nx * ny * nz * nw);
- if (index[6] >= nx * ny * nz * nw || index[6] < 0)
- fprintf(
- stderr,
- "index 6 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[6], i, j + 1, k + 1, l, nx * ny * nz * nw);
- if (index[7] >= nx * ny * nz * nw || index[7] < 0)
- fprintf(
- stderr,
- "index 7 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[7], i, j + 1, k + 1, l + 1, nx * ny * nz * nw);
- if (index[8] >= nx * ny * nz * nw || index[8] < 0)
- fprintf(
- stderr,
- "index 8 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[8], i + 1, j, k, l, nx * ny * nz * nw);
- if (index[9] >= nx * ny * nz * nw || index[9] < 0)
- fprintf(
- stderr,
- "index 9 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[9], i + 1, j, k, l + 1, nx * ny * nz * nw);
- if (index[10] >= nx * ny * nz * nw || index[10] < 0)
- fprintf(
- stderr,
- "index 10 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[10], i + 1, j, k + 1, l, nx * ny * nz * nw);
- if (index[11] >= nx * ny * nz * nw || index[11] < 0)
- fprintf(
- stderr,
- "index 11 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[11], i + 1, j, k + 1, l + 1, nx * ny * nz * nw);
- if (index[12] >= nx * ny * nz * nw || index[12] < 0)
- fprintf(
- stderr,
- "index 12 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[12], i + 1, j + 1, k, l, nx * ny * nz * nw);
- if (index[13] >= nx * ny * nz * nw || index[13] < 0)
- fprintf(
- stderr,
- "index 13 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[13], i + 1, j + 1, k, l + 1, nx * ny * nz * nw);
- if (index[14] >= nx * ny * nz * nw || index[14] < 0)
- fprintf(
- stderr,
- "index 14 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[14], i + 1, j + 1, k + 1, l, nx * ny * nz * nw);
- if (index[15] >= nx * ny * nz * nw || index[15] < 0)
- fprintf(
- stderr,
- "index 15 out of bounds %d, i,j,k,l, %d, %d, %d, %d, table size %d\n",
- index[15], i + 1, j + 1, k + 1, l + 1, nx * ny * nz * nw);
-#endif
-
- 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) {
- 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) {
- 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 2d EAGLE table in redshift and hydrogen
- * number density. Produces 1d table depending only
- * on log(temperature)
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_2d(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x,
- int n_x, int array_size, double *result_table, float *ub, float *lb) {
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- int index[2];
- for (int i = index_low; i < index_high; i++) {
- index[0] = row_major_index_2d(x_i, i, n_x,
- array_size);
- index[1] = row_major_index_2d(x_i + 1, i, n_x,
- array_size);
-
- result_table[i] = (1 - d_x) * table[index[0]] +
- d_x * table[index[1]];
- }
-}
-
-/*
- * @brief Interpolates 3d EAGLE table in redshift and hydrogen
- * number density. Produces 1d table depending only
- * on log(temperature)
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_3d(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x,
- int n_x, int y_i, float d_y, int n_y, int array_size, double *result_table, float *ub, float *lb) {
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- int index[4];
- for (int i = index_low; i < index_high; i++) {
- index[0] = row_major_index_3d(x_i, y_i, i, n_x,
- n_y, array_size);
- index[1] = row_major_index_3d(x_i, y_i + 1, i, n_x,
- n_y, array_size);
- index[2] = row_major_index_3d(x_i + 1, y_i, i, n_x,
- n_y, array_size);
- index[3] = row_major_index_3d(x_i + 1, y_i + 1, i, n_x,
- n_y, array_size);
-
- result_table[i] = (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]];
- }
-}
-
-/*
- * @brief Interpolates 4d EAGLE hydrogen and helium cooling
- * table in redshift, helium fraction and hydrogen number
- * density. Produces 1d table depending only on log(temperature)
- * and location of maximum cooling gradient with respect to
- * internal energy. NOTE: calculation of location of maximum
- * cooling gradient is still work in progress.
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param He_i Helium fraction index
- * @param d_He Helium fraction offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- * @param x_max_grad Pointer to location of maximum gradient
- * @param u_ini Internal energy at start of hydro timestep
- * @param dt Hydro timestep
- */
-__attribute__((always_inline)) INLINE static void
-construct_1d_table_from_4d_H_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 *table, int x_i, float d_x, int n_x,
- int y_i, float d_y, int n_y, int w_i, float d_w, int n_w, int array_size, double *result_table,
- double *x_max_grad, double u_ini, float dt, float *ub, float *lb) {
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- int index[8];
- for (int i = index_low; i < index_high; i++) {
- index[0] =
- row_major_index_4d(x_i, y_i, w_i, i, n_x,
- n_y, n_w, array_size);
- index[1] =
- row_major_index_4d(x_i, y_i, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[2] =
- row_major_index_4d(x_i, y_i + 1, w_i, i, n_x,
- n_y, n_w, array_size);
- index[3] =
- row_major_index_4d(x_i, y_i + 1, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[4] =
- row_major_index_4d(x_i + 1, y_i, w_i, i, n_x,
- n_y, n_w, array_size);
- index[5] =
- row_major_index_4d(x_i + 1, y_i, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[6] =
- row_major_index_4d(x_i + 1, y_i + 1, w_i, i, n_x,
- n_y, n_w, array_size);
- index[7] = row_major_index_4d(x_i + 1, y_i + 1, w_i + 1, i,
- n_x, n_y,
- n_w, array_size);
-
- result_table[i] = (1 - d_x) * (1 - d_y) * (1 - d_w) * table[index[0]] +
- (1 - d_x) * (1 - d_y) * d_w * table[index[1]] +
- (1 - d_x) * d_y * (1 - d_w) * table[index[2]] +
- (1 - d_x) * d_y * d_w * table[index[3]] +
- d_x * (1 - d_y) * (1 - d_w) * table[index[4]] +
- d_x * (1 - d_y) * d_w * table[index[5]] +
- d_x * d_y * (1 - d_w) * table[index[6]] +
- d_x * d_y * d_w * table[index[7]];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 1 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e "
- "%.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",
- i, d_x, d_y, d_w, table[index[0]], table[index[1]],
- table[index[2]], table[index[3]], table[index[4]], table[index[5]],
- table[index[6]], table[index[7]]);
-#endif
- }
-}
-
-/*
- * @brief Interpolates 4d EAGLE table in redshift,
- * helium fraction and hydrogen number density.
- * Produces 1d table depending only on log(temperature)
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param He_i Helium fraction index
- * @param d_He Helium fraction offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void construct_1d_table_from_4d(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x, int n_x,
- int y_i, float d_y, int n_y, int w_i, float d_w, int n_w, int array_size, double *result_table, float *ub, float *lb) {
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- int index[8];
- for (int i = index_low; i < index_high; i++) {
- index[0] =
- row_major_index_4d(x_i, y_i, w_i, i, n_x,
- n_y, n_w, array_size);
- index[1] =
- row_major_index_4d(x_i, y_i, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[2] =
- row_major_index_4d(x_i, y_i + 1, w_i, i, n_x,
- n_y, n_w, array_size);
- index[3] =
- row_major_index_4d(x_i, y_i + 1, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[4] =
- row_major_index_4d(x_i + 1, y_i, w_i, i, n_x,
- n_y, n_w, array_size);
- index[5] =
- row_major_index_4d(x_i + 1, y_i, w_i + 1, i, n_x,
- n_y, n_w, array_size);
- index[6] =
- row_major_index_4d(x_i + 1, y_i + 1, w_i, i, n_x,
- n_y, n_w, array_size);
- index[7] = row_major_index_4d(x_i + 1, y_i + 1, w_i + 1, i,
- n_x, n_y,
- n_w, array_size);
-
- result_table[i] = (1 - d_x) * (1 - d_y) * (1 - d_w) * table[index[0]] +
- (1 - d_x) * (1 - d_y) * d_w * table[index[1]] +
- (1 - d_x) * d_y * (1 - d_w) * table[index[2]] +
- (1 - d_x) * d_y * d_w * table[index[3]] +
- d_x * (1 - d_y) * (1 - d_w) * table[index[4]] +
- d_x * (1 - d_y) * d_w * table[index[5]] +
- d_x * d_y * (1 - d_w) * table[index[6]] +
- d_x * d_y * d_w * table[index[7]];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 2 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e "
- "%.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",
- i, d_x, d_y, d_w, table[index[0]], table[index[1]],
- table[index[2]], table[index[3]], table[index[4]], table[index[5]],
- table[index[6]], table[index[7]]);
-#endif
- }
-}
-
-/*
- * @brief calculates heating due to helium reionization
- *
- * @param z redshift
- * @param dz redshift offset
- */
-__attribute__((always_inline)) INLINE double eagle_helium_reionization_extraheat(double z, double dz, const struct cooling_function_data *restrict cooling) {
-
- double extra_heating = 0.0;
-
-/* dz is the change in redshift (start to finish) and hence *should* be < 0 */
-#ifdef SWIFT_DEBUG_CHECKS
- if (dz > 0) {
- error(" formulation of helium reionization expects dz<0, whereas you have "
- "dz=%e\n", dz);
- }
-#endif
-
- /* Helium reionization */
- if (cooling->he_reion == 1) {
- double he_reion_erg_pG = cooling->he_reion_ev_pH * eagle_ev_to_erg / eagle_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)) -
- erf((z - cooling->he_reion_z_center) /
- (pow(2., 0.5) * cooling->he_reion_z_sigma))) /
- 2.0;
- } else
- extra_heating = 0.0;
-
- return extra_heating;
-}
-
-/*
- * @brief Interpolates 4d EAGLE metal cooling tables in redshift,
- * helium fraction and hydrogen number density.
- * Produces 1d table depending only on log(temperature)
- * NOTE: Will need to be modified to incorporate particle to solar
- * electron abundance ratio.
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param He_i Helium fraction index
- * @param d_He Helium fraction offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void
-construct_1d_print_table_from_4d_elements(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x,
- int n_x, int y_i, float d_y, int n_y, int array_size, double *result_table, float *abundance_ratio, float *ub, float *lb) {
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- int index[4];
-
- for (int j = 0; j < cooling->N_Elements; j++) {
- for (int i = index_low; i < index_high; i++) {
- if (j == 0) result_table[i] = 0.0;
- index[0] = row_major_index_4d(x_i, y_i, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[1] = row_major_index_4d(x_i, y_i+1, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[2] = row_major_index_4d(x_i+1, y_i, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[3] = row_major_index_4d(x_i+1, y_i+1, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
-
- result_table[i+array_size*j] = ((1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]) * abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * (1 - d_y) * table[index[0]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]);
-#endif
- }
- }
-}
-
-/*
- * @brief Interpolates 4d EAGLE metal cooling tables in redshift,
- * helium fraction and hydrogen number density.
- * Produces 1d table depending only on log(temperature)
- * NOTE: Will need to be modified to incorporate particle to solar
- * electron abundance ratio.
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param He_i Helium fraction index
- * @param d_He Helium fraction offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void
-construct_1d_table_from_4d_elements(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x,
- int n_x, int y_i, float d_y, int n_y, int array_size, double *result_table, float *abundance_ratio, float *ub, float *lb) {
- int index[4];
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- for (int j = 0; j < cooling->N_Elements; j++) {
- for (int i = index_low; i < index_high; i++) {
- if (j == 0) result_table[i] = 0.0;
- index[0] = row_major_index_4d(x_i, y_i, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[1] = row_major_index_4d(x_i, y_i+1, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[2] = row_major_index_4d(x_i+1, y_i, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
- index[3] = row_major_index_4d(x_i+1, y_i+1, i, j, n_x,
- n_y, array_size, cooling->N_Elements);
-
- result_table[i] += ((1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]) * abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * (1 - d_y) * table[index[0]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]);
-#endif
- }
- }
-}
-
-/*
- * @brief Interpolates 3d EAGLE cooling tables in redshift,
- * helium fraction and hydrogen number density.
- * Produces 1d table depending only on log(temperature)
- * NOTE: Will need to be modified to incorporate particle to solar
- * electron abundance ratio.
- *
- * @param p Particle structure
- * @param cooling Cooling data structure
- * @param cosmo Cosmology structure
- * @param internal_const Physical constants structure
- * @param table Pointer to table we are interpolating
- * @param z_i Redshift index
- * @param d_z Redshift offset
- * @param He_i Helium fraction index
- * @param d_He Helium fraction offset
- * @param n_h_i Hydrogen number density index
- * @param d_n_h Hydrogen number density offset
- * @param result_table Pointer to 1d interpolated table
- */
-__attribute__((always_inline)) INLINE static void
-construct_1d_table_from_3d_elements(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const, float *table, int x_i, float d_x,
- int n_x, int array_size, double *result_table, float *abundance_ratio, float *ub, float *lb) {
- int index[2];
-
- const double log_10 = 2.302585092994;
- int index_low, index_high;
- float d_high, d_low;
-
- get_index_1d(cooling->Temp, cooling->N_Temp,(*ub)/log_10,&index_high,&d_high);
- get_index_1d(cooling->Temp, cooling->N_Temp,(*lb)/log_10,&index_low,&d_low);
- if (index_high < array_size-1) {
- index_high += 2;
- } else {
- index_high = array_size;
- }
- if (index_low > 0) index_low--;
-
- for (int j = 0; j < cooling->N_Elements; j++) {
- for (int i = index_low; i < index_high; i++) {
- if (j == 0) result_table[i] = 0.0;
- index[0] = row_major_index_3d(j, x_i, i,
- cooling->N_Elements, n_x,
- array_size);
- index[1] = row_major_index_3d(j, x_i + 1, i,
- cooling->N_Elements, n_x,
- array_size);
-
- result_table[i] += ((1 - d_x) * table[index[0]] +
- d_x * table[index[1]])*abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * table[index[0]],
- (1 - d_x) * table[index[0]] +
- d_x * table[index[1]]);
-#endif
- }
- }
-}
-
-/*
- * @brief Interpolates temperature from internal energy based on table
- *
- * @param temp Temperature
- * @param temperature_table Table of EAGLE temperature values
- * @param cooling Cooling data structure
- */
-__attribute__((always_inline)) INLINE static double
-eagle_convert_temp_to_u_1d_table(double temp, float *temperature_table,
- const struct cooling_function_data *restrict
- cooling) {
-
- int temp_i;
- float d_temp, u;
-
- get_index_1d(temperature_table, cooling->N_Temp, log10(temp), &temp_i,
- &d_temp);
-
- u = pow(10.0, interpol_1d(cooling->Therm, temp_i, d_temp));
-
- return u;
-}
-
-/*
- * @brief Interpolates temperature from internal energy based on table and
- * calculates the size of the internal energy cell for the specified
- * temperature.
- *
- * @param u Internal energy
- * @param delta_u Pointer to size of internal energy cell
- * @param temperature_table Table of EAGLE temperature values
- * @param cooling Cooling data structure
- */
-__attribute__((always_inline)) INLINE static 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 part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const) {
-
- int u_i;
- float d_u, logT;
- double upper, lower;
-
- get_index_1d(cooling->Therm, cooling->N_Temp, log_10_u, &u_i, &d_u);
-
- 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);
-
- *delta_u =
- pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
-
- return logT;
-}
-
-
-/*
- * @brief Interpolates temperature from internal energy based on table and
- * calculates the size of the internal energy cell for the specified
- * temperature.
- *
- * @param u Internal energy
- * @param delta_u Pointer to size of internal energy cell
- * @param temperature_table Table of EAGLE temperature values
- * @param cooling Cooling data structure
- */
-__attribute__((always_inline)) INLINE static double
-eagle_convert_u_to_temp_1d_table(
- double log_10_u, float *delta_u, double *temperature_table,
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const) {
-
- int u_i;
- float d_u, logT;//, T;
-
- get_index_1d(cooling->Therm, cooling->N_Temp, log_10_u, &u_i, &d_u);
-
- logT = interpol_1d_dbl(temperature_table, u_i, d_u);
-
- *delta_u =
- pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
-
- return logT;
-}
-
-/*
- * @brief Calculates cooling rate from multi-d tables
- *
- * @param u Internal energy
- * @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 element_lambda Pointer for printing contribution to cooling rate
- * from each of the metals.
- */
-__attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(
- double log_10_u, double *dlambda_du, 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,
- double *element_lambda,
- float *solar_ratio) {
- double T_gam, solar_electron_abundance, solar_electron_abundance1, solar_electron_abundance2, elem_cool1, elem_cool2;
- 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, temp_lambda1, temp_lambda2;
- float du;
- double h_plus_he_electron_abundance;
- double h_plus_he_electron_abundance1;
- double h_plus_he_electron_abundance2;
-
- int i;
- double temp;
- int temp_i;
- float d_temp;
-
- *dlambda_du = 0.0;
-
- temp = eagle_convert_u_to_temp(log_10_u, &du, z_index, n_h_i, He_i,
- dz, d_n_h, d_He, p,
- cooling, cosmo, internal_const);
-
- get_index_1d(cooling->Temp, cooling->N_Temp, temp, &temp_i, &d_temp);
-
- /* ------------------ */
- /* Metal-free cooling */
- /* ------------------ */
-
- /* redshift tables */
- temp_lambda = interpol_4d(cooling->table.element_cooling.H_plus_He_heating,
- z_index, n_h_i, He_i, temp_i, dz, d_n_h, d_He,
- d_temp, cooling->N_Redshifts, cooling->N_nH,
- cooling->N_He, cooling->N_Temp,&temp_lambda2,&temp_lambda1);
- h_plus_he_electron_abundance = interpol_4d(
- cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- n_h_i, He_i, temp_i, dz, d_n_h, d_He, d_temp, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp,&h_plus_he_electron_abundance2,&h_plus_he_electron_abundance1);
- cooling_rate += temp_lambda;
- *dlambda_du += (temp_lambda2 - temp_lambda1)/du;
- if (element_lambda != NULL) element_lambda[0] = temp_lambda;
-
- /* ------------------ */
- /* Compton cooling */
- /* ------------------ */
-
- 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 */
+double eagle_metal_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 *,
+ double *,
+ float *);
- 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, find the abundance and multiply it
- by the interpolated heating-cooling */
-
- /* redshift tables */
- solar_electron_abundance =
- interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,
- n_h_i, temp_i, dz, d_n_h, d_temp, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp,&solar_electron_abundance2,&solar_electron_abundance1);
-
- for (i = 0; i < cooling->N_Elements; i++) {
- temp_lambda =
- interpol_4d(cooling->table.element_cooling.metal_heating, z_index,
- n_h_i, temp_i, i, dz, d_n_h, d_temp, 0.0, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp, cooling->N_Elements,&elem_cool2,&elem_cool1) *
- (h_plus_he_electron_abundance / solar_electron_abundance)*
- solar_ratio[i+2];
- cooling_rate += temp_lambda;
- *dlambda_du +=
- (elem_cool2 * h_plus_he_electron_abundance2 /
- solar_electron_abundance2 -
- elem_cool1 * h_plus_he_electron_abundance1 /
- solar_electron_abundance1) / du * solar_ratio[i+2];
- if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
- }
-
- return cooling_rate;
-}
-
-/*
- * @brief Calculates cooling rate from 1d tables
- *
- * @param u Internal energy
- * @param dlambda_du Pointer to gradient of 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 element_lambda Pointer for printing contribution to cooling rate
- * from each of the metals.
- */
-__attribute__((always_inline)) INLINE static double
+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, 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,
- double *element_lambda,
- float *solar_ratio) {
- 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;
- double temp;
- int temp_i;
- float d_temp;
-
- *dlambda_du = 0.0;
-
- temp = eagle_convert_u_to_temp_1d_table(log_10_u, &du, temp_table, p, cooling, cosmo,
- internal_const);
-
- get_index_1d(cooling->Temp, cooling->N_Temp, temp, &temp_i, &d_temp);
- //printf("Eagle cooling.h 1d log_10_u, temp, temp_i, d_temp %.5e %.5e %d %.5e\n", log_10_u, temp, temp_i, d_temp);
-
- /* ------------------ */
- /* Metal-free cooling */
- /* ------------------ */
-
- /* redshift tables */
- 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 */
- /* ------------------ */
-
- 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, find the abundance and multiply it
- by the interpolated heating-cooling */
-
- /* redshift 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); // hydrogen and helium aren't included here.
- 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 previously used to calculate cooling rate and
- * dLambda_du
- * Can be used to check whether calculation of dLambda_du is done correctly.
- * Should be removed when finished
- *
- * @param u Internal energy
- * @param dlambda_du Pointer to gradient of 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
- */
-__attribute__((always_inline)) INLINE static double eagle_cooling_rate(
- double logu, double *dLambdaNet_du, 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;
- const double log_10 = 2.30258509299404;
-
- lambda_net = eagle_metal_cooling_rate(
- logu/log_10, dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
-
- return lambda_net;
-}
-
-/*
- * @brief Wrapper function previously used to calculate cooling rate and
- * dLambda_du
- * Can be used to check whether calculation of dLambda_du is done correctly.
- * Should be removed when finished
- *
- * @param u Internal energy
- * @param dlambda_du Pointer to gradient of 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
- */
-__attribute__((always_inline)) INLINE static 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;
- const double log_10 = 2.30258509299404;
-
- lambda_net = eagle_metal_cooling_rate_1d_table(
- logu/log_10, dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_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, element_lambda, abundance_ratio);
-
- return lambda_net;
-}
-
-/*
- * @brief Calculates cooling rate from 1d tables, used to print cooling rates
- * to files for checking against Wiersma et al 2009 idl scripts.
- *
- * @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
- */
-__attribute__((always_inline)) INLINE static double
+ 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 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, lambda_net = 0.0, dLambdaNet_du;
- element_lambda = malloc((cooling->N_Elements + 2) * sizeof(double));
- double u = hydro_get_physical_internal_energy(p, cosmo) *
- cooling->internal_energy_scale;
-
- 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(
- log10(u), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
- 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 Calculates cooling rate from 1d tables, used to print cooling rates
- * to files for checking against Wiersma et al 2009 idl scripts.
- *
- * @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
- */
-__attribute__((always_inline)) INLINE static 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_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, 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, 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, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
- 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 for when Newton Raphson method fails to
- * converge
- *
- * @param x_init Initial guess for log(internal energy)
- * @param u_ini Internal energy at beginning of hydro step
- * @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 cosmo Cosmology structure
- * @param cooling Cooling data structure
- * @param phys_const Physical constants structure
- * @param dt Hydro timestep
- */
-__attribute__((always_inline)) INLINE static float bisection_iter(
- float x_init, double u_ini, 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, float He_reion_heat,
- struct part *restrict p, const struct cosmology *restrict cosmo,
- const struct cooling_function_data *restrict cooling,
- const struct phys_const *restrict phys_const,
- float *abundance_ratio, float dt, float *ub, float *lb) {
- double x_a, x_b, x_next, f_a, f_b, f_next, dLambdaNet_du;
- int i = 0;
- const float bisection_tolerance = 1.0e-8;
-
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-
- /* convert Hydrogen mass fraction in Hydrogen number density */
- double inn_h =
- chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
- cooling->number_density_scale;
-
- /* 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);
-
- /* set helium and hydrogen reheating term */
-
- // Add bracketing?
- x_a = *ub;
- x_b = *lb;
- f_a = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate_1d_table(
- x_a, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- f_b = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate_1d_table(
- x_b, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
-
-#if SWIFT_DEBUG_CHECKS
- assert(f_a * f_b < 0);
-#endif
-
- i = 0;
- do {
- x_next = 0.5 * (x_a + x_b);
- f_next = (He_reion_heat / (dt * ratefact)) + eagle_cooling_rate_1d_table(
- x_next, &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_table,
- element_electron_abundance_table, temp_table, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- if (f_a * f_next < 0) {
- x_b = x_next;
- f_b = f_next;
- } else {
- x_a = x_next;
- f_a = f_next;
- }
- i++;
- } while (fabs(x_a - x_b) > bisection_tolerance && i < 50*eagle_max_iterations);
-
- if (i >= 50*eagle_max_iterations) n_eagle_cooling_rate_calls_4++;
-
- return x_b;
-}
-
-
-/*
- * @brief Newton Raphson integration scheme to calculate particle cooling over
- * hydro timestep
- *
- * @param x_init Initial guess for log(internal energy)
- * @param u_ini Internal energy at beginning of hydro step
- * @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 cosmo Cosmology structure
- * @param cooling Cooling data structure
- * @param phys_const Physical constants structure
- * @param dt Hydro timestep
- * @param printflag Flag to identify if scheme failed to converge
- * @param lb lower bound to be used by bisection scheme
- * @param ub upper bound to be used by bisection scheme
- */
-__attribute__((always_inline)) INLINE static float newton_iter(
- float x_init, double u_ini, int z_index,
- float dz, int n_h_i, float d_n_h, int He_i, float d_He,
- float He_reion_heat, struct part *restrict p,
- const struct cosmology *restrict cosmo,
- const struct cooling_function_data *restrict cooling,
- const struct phys_const *restrict phys_const,
- float *abundance_ratio, float dt, int *printflag,
- float *lb, float *ub) {
- /* this routine does the iteration scheme, call one and it iterates to
- * convergence */
-
- const float newton_tolerance = 1.0e-2;
- double x, x_old;
- double dLambdaNet_du, LambdaNet;
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
-
- /* convert Hydrogen mass fraction in Hydrogen number density */
- double inn_h =
- chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
- cooling->number_density_scale;
-
- /* 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);
-
- x_old = x_init;
- x = x_old;
- int i = 0;
- float relax = 1.0;
- float log_table_bound_high = 43.749, log_table_bound_low = 20.723; // Corresponds to u = 10^19, 10^9
-
- do /* iterate to convergence */
- {
- if (relax == 1.0) x_old = x;
- LambdaNet = (He_reion_heat / (dt * ratefact)) +
- eagle_cooling_rate(
- x_old, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
- if (LambdaNet < 0) {
- *ub = fmin(*ub,x_old);
- } else if (LambdaNet > 0) {
- *lb = fmax(*lb,x_old);
- }
-
- n_eagle_cooling_rate_calls_1++;
-
- // Newton iteration.
- x = x_old -
- relax * (1.0 - u_ini * exp(-x_old) -
- LambdaNet * ratefact * dt * exp(-x_old)) /
- (1.0 - dLambdaNet_du * ratefact * dt);
-
- // check whether iterations go out of bounds of table,
- // if out of bounds, try again with smaller newton step
- if (x > log_table_bound_high) {
- x = (log_table_bound_high + x_old)/2.0;
- } else if (x < log_table_bound_low) {
- x = (log_table_bound_low + x_old)/2.0;
- }
-
- i++;
- } while (fabs(x - x_old) > newton_tolerance && i < eagle_max_iterations);
- if (i >= eagle_max_iterations) {
- n_eagle_cooling_rate_calls_3++;
- // failed to converge the first time
- if (*printflag == 0) *printflag = 1;
- }
-
- return x;
-}
-
-/*
- * @brief Calculate ratio of particle element abundances
- * to solar abundance
- *
- * @param p Pointer to particle data
- * @param cooling Pointer to cooling data
- */
-__attribute__((always_inline)) INLINE static void abundance_ratio_to_solar(
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- float *ratio_solar) {
- for(enum chemistry_element elem = chemistry_element_H; elem < chemistry_element_count; elem++){
- if (elem == chemistry_element_Fe) {
- ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem]/
- cooling->SolarAbundances[elem+2]; // NOTE: solar abundances have iron last with calcium and sulphur directly before, hence +2
- } else {
- ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem]/
- cooling->SolarAbundances[elem];
- }
- }
- ratio_solar[chemistry_element_count] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si]*
- cooling->sulphur_over_silicon_ratio/
- cooling->SolarAbundances[chemistry_element_count-1]; // NOTE: solar abundances have iron last with calcium and sulphur directly before, hence -1
- ratio_solar[chemistry_element_count+1] = p->chemistry_data.metal_mass_fraction[chemistry_element_Si]*
- cooling->calcium_over_silicon_ratio/
- cooling->SolarAbundances[chemistry_element_count]; // NOTE: solar abundances have iron last with calcium and sulphur directly before
-}
-
-/**
- * @brief construct all 1d tables
- *
- * @param
- */
-__attribute__((always_inline)) INLINE void static construct_1d_tables(
- int z_index, float dz, int n_h_i, float d_n_h,
- int He_i, float d_He,
- const struct phys_const *restrict phys_const,
- const struct unit_system *restrict us,
- const struct cosmology *restrict cosmo,
- const struct cooling_function_data *restrict cooling,
- const struct part *restrict p,
- float *abundance_ratio,
- double *H_plus_He_heat_table,
- double *H_plus_He_electron_abundance_table,
- double *element_cooling_table,
- double *element_cooling_print_table,
- double *element_electron_abundance_table,
- double *temp_table,
- float *ub, float *lb) {
-
- construct_1d_table_from_4d(p, cooling, cosmo, phys_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);
-
- // element cooling
- construct_1d_table_from_4d(
- p, cooling, cosmo, phys_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, phys_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, phys_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, phys_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_print_table, abundance_ratio, ub, lb);
- construct_1d_table_from_3d(
- p, cooling, cosmo, phys_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);
-
- // photodissociation cooling
- //construct_1d_table_from_3d(
- // p, cooling, cosmo, phys_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, phys_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, phys_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_table_from_2d(
- // p, cooling, cosmo, phys_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);
-}
-
-/**
- * @brief Apply the cooling function to a particle.
- *
- * @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.
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- * @param dt The time-step of this particle.
- */
-__attribute__((always_inline)) INLINE static void cooling_cool_part(
- const struct phys_const *restrict phys_const,
- const struct unit_system *restrict us,
- const struct cosmology *restrict cosmo,
- const struct cooling_function_data *restrict cooling,
- struct part *restrict p, struct xpart *restrict xp, float dt) {
-
- const float exp_factor = 0.05;
- const double log_10_e = 0.43429448190325182765;
- double u_old = (hydro_get_physical_internal_energy(p, cosmo)
- + hydro_get_internal_energy_dt(p) * dt) *
- cooling->internal_energy_scale;
- float dz;
- int z_index;
- get_redshift_index(cosmo->z, &z_index, &dz, cooling);
-
- float XH, HeFrac;
- double inn_h;
-
- double ratefact, u, LambdaNet, LambdaTune = 0, dLambdaNet_du;
-
- dt *= units_cgs_conversion_factor(us, UNIT_CONV_TIME);
-
- u = u_old;
- double u_ini = u_old;//, u_temp;
-
- float abundance_ratio[chemistry_element_count + 2];
- abundance_ratio_to_solar(p, cooling, abundance_ratio);
-
- 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]);
-
- /* convert Hydrogen mass fraction in Hydrogen number density */
- inn_h =
- chemistry_get_number_density(p, cosmo, chemistry_element_H, phys_const) *
- cooling->number_density_scale;
- /* 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);
- /* set helium and hydrogen reheating term */
- if (cooling->he_reion == 1) LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
-
- int He_i, n_h_i;
- float d_He, d_n_h;
- 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);
-
-
- LambdaNet = LambdaTune/(dt*ratefact) + eagle_cooling_rate(
- log(u_ini), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
-
- if (fabs(ratefact * LambdaNet * dt) < exp_factor * u_old) {
- /* cooling rate is small, take the explicit solution */
- u = u_old + ratefact * LambdaNet * dt;
- } else {
- n_eagle_cooling_rate_calls_2++;
- float x, lower_bound = cooling->Therm[0]/log_10_e, upper_bound = cooling->Therm[cooling->N_Temp-1]/log_10_e;
-
- double u_eq = u_ini;
-
- if (dt > 0) {
- double log_u_temp = log(u_eq);
-
- int printflag = 0;
- x = newton_iter(log_u_temp, u_ini, z_index, dz,
- n_h_i, d_n_h, He_i, d_He, LambdaTune, p, cosmo, cooling, phys_const,
- abundance_ratio, dt, &printflag, &lower_bound, &upper_bound);
- if (printflag == 1) {
- // newton method didn't work, so revert to bisection
-
- // construct 1d table of cooling rates wrt temperature
- double temp_table[176]; // WARNING sort out how it is declared/allocated
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double element_cooling_table[176];
- double element_cooling_print_table[9*176];
- double element_electron_abundance_table[176];
- construct_1d_tables(z_index, dz, n_h_i, d_n_h, He_i, d_He,
- phys_const, us, cosmo, cooling, p,
- abundance_ratio,
- H_plus_He_heat_table,
- H_plus_He_electron_abundance_table,
- element_cooling_table,
- element_cooling_print_table,
- element_electron_abundance_table,
- temp_table,
- &upper_bound,
- &lower_bound);
- // perform bisection scheme
- x = bisection_iter(log_u_temp,u_ini,H_plus_He_heat_table,H_plus_He_electron_abundance_table,element_cooling_print_table,element_electron_abundance_table,temp_table,z_index,dz,n_h_i,d_n_h,He_i,d_He,LambdaTune,p,cosmo,cooling,phys_const,abundance_ratio,dt,&upper_bound,&lower_bound);
-
- printflag = 2;
- }
- u = exp(x);
- }
- }
-
- const float hydro_du_dt = hydro_get_internal_energy_dt(p);
-
-
- // calculate du/dt
- float cooling_du_dt = 0.0;
- if (dt > 0) {
- cooling_du_dt = (u - u_ini) / dt / cooling->power_scale;
- }
-
- /* Update the internal energy time derivative */
- hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
-
- /* Store the radiated energy */
- xp->cooling_data.radiated_energy +=
- -hydro_get_mass(p) * cooling_du_dt *
- (dt / units_cgs_conversion_factor(us, UNIT_CONV_TIME));
-}
-
-/**
- * @brief Writes the current model of SPH to the file
- * @param h_grpsph The HDF5 group in which to write
- */
-__attribute__((always_inline)) INLINE static void cooling_write_flavour(
- hid_t h_grpsph) {
-
- io_write_attribute_s(h_grpsph, "Cooling Model", "EAGLE");
-}
-
-/**
- * @brief Computes the cooling time-step.
- *
- * @param cooling The #cooling_function_data used in the run.
- * @param phys_const The physical constants in internal units.
- * @param us The internal system of units.
- * @param cosmo The current cosmological model.
- * @param p Pointer to the particle data.
- */
-__attribute__((always_inline)) INLINE static float cooling_timestep(
- const struct cooling_function_data *restrict cooling,
- const struct phys_const *restrict phys_const,
- const struct cosmology *restrict cosmo,
- const struct unit_system *restrict us, const struct part *restrict p, const struct xpart *restrict xp) {
-
- /* Remember to update when using an implicit integrator!!!*/
- // const float cooling_rate = cooling->cooling_rate;
- // const float internal_energy = hydro_get_comoving_internal_energy(p);
- // return cooling->cooling_tstep_mult * internal_energy / fabsf(cooling_rate);
-
- return FLT_MAX;
-}
-
-/**
- * @brief Sets the cooling properties of the (x-)particles to a valid start
- * state.
- *
- * @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 properties of the cooling function.
- * @param p Pointer to the particle data.
- * @param xp Pointer to the extended particle data.
- */
-__attribute__((always_inline)) INLINE static void cooling_first_init_part(
- const struct phys_const* restrict phys_const,
- const struct unit_system* restrict us,
- const struct cosmology* restrict cosmo,
- const struct cooling_function_data* restrict cooling,
- const struct part* restrict p, struct xpart* restrict xp) {
-
- xp->cooling_data.radiated_energy = 0.f;
-}
-
-/**
- * @brief Returns the total radiated energy by this particle.
- *
- * @param xp The extended particle data
- */
-__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
- const struct xpart *restrict xp) {
-
- return xp->cooling_data.radiated_energy;
-}
-
-/**
- * @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 static 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");
- // MATTHIEU: 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);
- }
-}
-
-/**
- * @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 static 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.
- *
- * @param parameter_file The parsed parameter file.
- * @param us The current internal system of units.
- * @param phys_const The physical constants in internal units.
- * @param cooling The cooling properties to initialize
- */
-static INLINE void cooling_init_backend(
- const struct swift_params *parameter_file, const struct unit_system *us,
- const struct phys_const *phys_const,
- struct cooling_function_data *cooling) {
-
- char fname[200];
-
- parser_get_param_string(parameter_file, "EagleCooling:filename",
- cooling->cooling_table_path);
- cooling->reionisation_redshift = parser_get_param_float(
- parameter_file, "EagleCooling:reionisation_redshift");
- cooling->calcium_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:CalciumOverSilicon");
- cooling->sulphur_over_silicon_ratio = parser_get_param_float(
- parameter_file, "EAGLEChemistry:SulphurOverSilicon");
- cooling->he_reion = 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");
- cooling->he_reion_z_sigma = parser_get_param_float(
- parameter_file, "EagleCooling:he_reion_z_sigma");
- cooling->he_reion_ev_pH = parser_get_param_float(
- parameter_file, "EagleCooling:he_reion_ev_pH");
-
- GetCoolingRedshifts(cooling);
- sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
- ReadCoolingHeader(fname, cooling);
- MakeNamePointers(cooling);
- cooling->table = eagle_readtable(cooling->cooling_table_path, cooling);
- printf("Eagle cooling.h read table \n");
-
- cooling->ElementAbundance_SOLARM1 =
- malloc(cooling->N_SolarAbundances * sizeof(float));
- cooling->solar_abundances = malloc(cooling->N_Elements * sizeof(double));
-
- cooling->delta_u = cooling->Therm[1] - cooling->Therm[0];
-
- cooling->internal_energy_scale =
- units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) /
- units_cgs_conversion_factor(us, UNIT_CONV_MASS);
- cooling->number_density_scale =
- units_cgs_conversion_factor(us, UNIT_CONV_DENSITY) /
- units_cgs_conversion_factor(us, UNIT_CONV_MASS);
- cooling->power_scale = units_cgs_conversion_factor(us, UNIT_CONV_POWER) /
- units_cgs_conversion_factor(us, UNIT_CONV_MASS);
- cooling->temperature_scale =
- units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
-}
-
-/**
- * @brief Prints the properties of the cooling model to stdout.
- *
- * @param cooling The properties of the cooling function.
- */
-static INLINE void cooling_print_backend(
- const struct cooling_function_data *cooling) {
-
- message("Cooling function is 'EAGLE'.");
-}
+ 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,
+ const struct phys_const *restrict,
+ float *, float);
+
+float newton_iter(
+ float, double, int,
+ float, int, float, int, float,
+ float, struct part *restrict,
+ const struct cosmology *restrict,
+ const struct cooling_function_data *restrict,
+ const struct phys_const *restrict,
+ float *, float, int *);
+
+void abundance_ratio_to_solar(
+ const struct part *restrict,
+ const struct cooling_function_data *restrict,
+ float *);
+
+void cooling_cool_part(
+ const struct phys_const *restrict,
+ const struct unit_system *restrict,
+ const struct cosmology *restrict,
+ const struct cooling_function_data *restrict,
+ struct part *restrict, struct xpart *restrict, float);
+
+void cooling_write_flavour(
+ hid_t);
+
+float cooling_timestep(
+ const struct cooling_function_data *restrict,
+ const struct phys_const *restrict,
+ const struct cosmology *restrict,
+ const struct unit_system *restrict, const struct part *restrict, const struct xpart *restrict);
+
+void cooling_first_init_part(
+ const struct phys_const* restrict,
+ const struct unit_system* restrict,
+ const struct cosmology* restrict,
+ const struct cooling_function_data* restrict,
+ const struct part* restrict, struct xpart* 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 *);
+
+void cooling_print_backend(const struct cooling_function_data *);
#endif /* SWIFT_COOLING_EAGLE_H */
diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h
index 97be7b11be56d62a7c6f1e58d5bf78fc157a4a17..4e903d580365f6b93f8da2132f26eedfd64ab87e 100644
--- a/src/cooling/EAGLE/cooling_struct.h
+++ b/src/cooling/EAGLE/cooling_struct.h
@@ -29,6 +29,8 @@
#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
diff --git a/src/cooling/EAGLE/eagle_cool_tables.c b/src/cooling/EAGLE/eagle_cool_tables.c
index 4aa8d8496095ebd703d4ecf5057c7037223aca25..58e423a58ef41e60ad671729adf88087d389b3b2 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.c
+++ b/src/cooling/EAGLE/eagle_cool_tables.c
@@ -31,10 +31,23 @@
#include <math.h>
#include <stdlib.h>
#include <string.h>
-#include "cooling_struct.h"
#include "eagle_cool_tables.h"
-#include "error.h"
+#include "cooling_struct.h"
#include "interpolate.h"
+#include "error.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
@@ -72,7 +85,8 @@ void GetCoolingRedshifts(struct cooling_function_data *cooling) {
* @param fname Filepath for cooling table from which to read header
* @param cooling Cooling data structure
*/
-void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
+void ReadCoolingHeader(char *fname,
+ struct cooling_function_data *cooling) {
int i;
hid_t tempfile_id, dataset, datatype;
@@ -121,14 +135,10 @@ 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 * 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));
+ cooling->ElementNames =
+ malloc(cooling->N_Elements * eagle_element_name_length * sizeof(char));
+ cooling->SolarAbundanceNames = malloc(
+ cooling->N_SolarAbundances * eagle_element_name_length * sizeof(char));
// read in values for each of the arrays
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins", H5P_DEFAULT);
@@ -172,12 +182,11 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
H5Tclose(datatype);
for (i = 0; i < cooling->N_Elements; i++)
- strcpy(cooling->ElementNames[i],element_names[i]);
-
+ cooling->ElementNames[i] = mystrdup(element_names[i]);
char solar_abund_names[cooling->N_SolarAbundances][eagle_element_name_length];
- // read in assumed solar abundances used in constructing the tables,
+ // read in assumed solar abundances used in constructing the tables,
// and corresponding names
datatype = H5Tcopy(H5T_C_S1);
status = H5Tset_size(datatype, string_length);
@@ -188,11 +197,11 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
H5Tclose(datatype);
for (i = 0; i < cooling->N_SolarAbundances; i++)
- strcpy(cooling->SolarAbundanceNames[i],element_names[i]);
+ cooling->SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
status = H5Fclose(tempfile_id);
- // Convert to temperature, density and internal energy arrays to log10
+ // 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]);
@@ -202,15 +211,14 @@ void ReadCoolingHeader(char *fname, struct cooling_function_data *cooling) {
}
/*
- * @brief Get the table of cooling rates for photoionized cooling (before
- * redshift ~9) Reads in table of cooling rates and electron abundances due to
- * metals (depending on temperature, hydrogen number density), cooling rates and
- * electron abundances due to hydrogen and helium (depending on temperature,
- * hydrogen number density and helium fraction), and temperatures (depending on
- * internal energy, hydrogen number density and helium fraction; note: this is
- * distinct from table of temperatures read in ReadCoolingHeader, as that table
- * is used to index the cooling, electron abundance tables, whereas this one is
- * used to obtain temperature of particle)
+ * @brief Get the table of cooling rates for photoionized cooling (before redshift ~9)
+ * Reads in table of cooling rates and electron abundances due to metals
+ * (depending on temperature, hydrogen number density), cooling rates and electron
+ * abundances due to hydrogen and helium (depending on temperature, hydrogen number
+ * density and helium fraction), and temperatures (depending on internal energy,
+ * hydrogen number density and helium fraction; note: this is distinct from table of
+ * temperatures read in ReadCoolingHeader, as that table is used to index the cooling,
+ * electron abundance tables, whereas this one is used to obtain temperature of particle)
*
* @param cooling_table_path Directory containing cooling table
* @param cooling Cooling data structure
@@ -235,7 +243,7 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
float *he_net_cooling_rate;
float *he_electron_abundance;
- // allocate temporary arrays (needed to change order of dimensions
+ // allocate temporary arrays (needed to change order of dimensions
// of arrays)
net_cooling_rate =
(float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
@@ -276,16 +284,16 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
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!!!
+ 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];
}
}
}
- // read in cooling rates due to hydrogen and helium, H + He electron
- // abundances, temperatures
+ // 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,
@@ -333,7 +341,8 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
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);
+ row_major_index_3d(0, j, i, 1, cooling->N_nH,
+ cooling->N_Temp);
cooling_table.electron_abundance[cooling_index] =
electron_abundance[table_index];
}
@@ -347,21 +356,20 @@ struct cooling_tables_redshift_invariant get_no_compt_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- message("done reading in no compton table");
+ printf("eagle_cool_tables.h done reading in no compton table\n");
return cooling_table;
}
/*
- * @brief Get the table of cooling rates for collisional cooling
+ * @brief Get the table of cooling rates for collisional cooling
* Reads in table of cooling rates and electron abundances due to metals
- * (depending on temperature, hydrogen number density), cooling rates and
- * electron abundances due to hydrogen and helium (depending on temperature,
- * hydrogen number density and helium fraction), and temperatures (depending on
- * internal energy, hydrogen number density and helium fraction; note: this is
- * distinct from table of temperatures read in ReadCoolingHeader, as that table
- * is used to index the cooling, electron abundance tables, whereas this one is
- * used to obtain temperature of particle)
+ * (depending on temperature, hydrogen number density), cooling rates and electron
+ * abundances due to hydrogen and helium (depending on temperature, hydrogen number
+ * density and helium fraction), and temperatures (depending on internal energy,
+ * hydrogen number density and helium fraction; note: this is distinct from table of
+ * temperatures read in ReadCoolingHeader, as that table is used to index the cooling,
+ * electron abundance tables, whereas this one is used to obtain temperature of particle)
*
* @param cooling_table_path Directory containing cooling table
* @param cooling Cooling data structure
@@ -385,7 +393,7 @@ struct cooling_tables_redshift_invariant get_collisional_table(
float *he_net_cooling_rate;
float *he_electron_abundance;
- // allocate temporary arrays (needed to change order of dimensions
+ // allocate temporary arrays (needed to change order of dimensions
// of arrays)
net_cooling_rate = (float *)malloc(cooling->N_Temp * sizeof(float));
electron_abundance = (float *)malloc(cooling->N_Temp * sizeof(float));
@@ -428,8 +436,8 @@ struct cooling_tables_redshift_invariant get_collisional_table(
}
}
- // read in cooling rates due to hydrogen and helium, H + He electron
- // abundances, temperatures
+ // 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,
@@ -482,20 +490,19 @@ struct cooling_tables_redshift_invariant get_collisional_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- message("done reading in collisional table");
+ printf("eagle_cool_tables.h done reading in collisional table\n");
return cooling_table;
}
/*
- * @brief Get the table of cooling rates for photodissociation cooling
+ * @brief Get the table of cooling rates for photodissociation cooling
* Reads in table of cooling rates and electron abundances due to metals
- * (depending on temperature, hydrogen number density), cooling rates and
- * electron abundances due to hydrogen and helium (depending on temperature,
- * hydrogen number density and helium fraction), and temperatures (depending on
- * internal energy, hydrogen number density and helium fraction; note: this is
- * distinct from table of temperatures read in ReadCoolingHeader, as that table
- * is used to index the cooling, electron abundance tables, whereas this one is
- * used to obtain temperature of particle)
+ * (depending on temperature, hydrogen number density), cooling rates and electron
+ * abundances due to hydrogen and helium (depending on temperature, hydrogen number
+ * density and helium fraction), and temperatures (depending on internal energy,
+ * hydrogen number density and helium fraction; note: this is distinct from table of
+ * temperatures read in ReadCoolingHeader, as that table is used to index the cooling,
+ * electron abundance tables, whereas this one is used to obtain temperature of particle)
*
* @param cooling_table_path Directory containing cooling table
* @param cooling Cooling data structure
@@ -519,7 +526,7 @@ struct cooling_tables_redshift_invariant get_photodis_table(
float *he_net_cooling_rate;
float *he_electron_abundance;
- // allocate temporary arrays (needed to change order of dimensions
+ // allocate temporary arrays (needed to change order of dimensions
// of arrays)
net_cooling_rate =
(float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
@@ -560,16 +567,16 @@ struct cooling_tables_redshift_invariant get_photodis_table(
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!!!
+ 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];
}
}
}
- // read in cooling rates due to hydrogen and helium, H + He electron
- // abundances, temperatures
+ // 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,
@@ -631,7 +638,7 @@ struct cooling_tables_redshift_invariant get_photodis_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- message("done reading in photodissociation table");
+ printf("eagle_cool_tables.h done reading in photodissociation table\n");
return cooling_table;
}
@@ -661,7 +668,7 @@ struct cooling_tables get_cooling_table(
float *he_net_cooling_rate;
float *he_electron_abundance;
- // allocate temporary arrays (needed to change order of dimensions
+ // allocate temporary arrays (needed to change order of dimensions
// of arrays)
net_cooling_rate =
(float *)malloc(cooling->N_Temp * cooling->N_nH * sizeof(float));
@@ -713,16 +720,16 @@ struct cooling_tables get_cooling_table(
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);
+ 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];
}
}
}
- // read in cooling rates due to hydrogen and helium, H + He electron
- // abundances, temperatures
+ // 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,
@@ -785,7 +792,7 @@ struct cooling_tables get_cooling_table(
free(he_net_cooling_rate);
free(he_electron_abundance);
- message("done reading in general cooling table");
+ printf("eagle_cool_tables.h done reading in general cooling table\n");
return cooling_table;
}
@@ -802,7 +809,8 @@ struct eagle_cooling_table eagle_readtable(
struct eagle_cooling_table table;
- table.no_compton_cooling = get_no_compt_table(cooling_table_path, cooling);
+ 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 =
@@ -812,4 +820,180 @@ 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/interpolate.c b/src/cooling/EAGLE/interpolate.c
index 1248f6c2881ea1b891822bf86fb5d2a578e6e04d..5bba67cbcb9cfe242efff42ce3026b97f0ae2cee 100644
--- a/src/cooling/EAGLE/interpolate.c
+++ b/src/cooling/EAGLE/interpolate.c
@@ -57,12 +57,7 @@ const float EPS = 1.e-4;
*/
__attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
int nx, int ny) {
- int index = i * ny + j;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
-#endif
- return index;
+ return i * ny + j;
}
/**
@@ -75,13 +70,7 @@ __attribute__((always_inline)) INLINE int row_major_index_2d(int i, int j,
__attribute__((always_inline)) INLINE int row_major_index_3d(int i, int j,
int k, int nx,
int ny, int nz) {
- int index = i * ny * nz + j * nz + k;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
- assert(k < nz);
-#endif
- return index;
+ return i * ny * nz + j * nz + k;
}
/**
@@ -95,14 +84,7 @@ __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) {
- int index = i * ny * nz * nw + j * nz * nw + k * nw + l;
-#ifdef SWIFT_DEBUG_CHECKS
- assert(i < nx);
- assert(j < ny);
- assert(k < nz);
- assert(l < nw);
-#endif
- return index;
+ return i * ny * nz * nw + j * nz * nw + k * nw + l;
}
/*
@@ -143,7 +125,7 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table,
*
* @param z Redshift whose position within the redshift array we are interested
* in
- * @param z_index i Pointer to the index whose corresponding redshift
+ * @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
@@ -229,6 +211,10 @@ __attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table,
* @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,
@@ -268,6 +254,10 @@ __attribute__((always_inline)) INLINE float interpol_2d(float *table, int i,
* @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 double interpol_2d_dbl(
double *table, int i, int j, double dx, double dy, int nx, int ny,
@@ -305,6 +295,10 @@ __attribute__((always_inline)) INLINE double interpol_2d_dbl(
* @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,
@@ -370,6 +364,14 @@ __attribute__((always_inline)) INLINE float interpol_3d(float *table, int i,
* @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,
@@ -428,6 +430,7 @@ __attribute__((always_inline)) INLINE float interpol_4d(
dx * dy * dz * (1 - dw) * table14 +
dx * dy * dz * dw * table15;
if (nw == 9) {
+ // interpolating metal species
dz = 1.0;
} else {
dw = 1.0;
@@ -449,6 +452,7 @@ __attribute__((always_inline)) INLINE float interpol_4d(
dx * dy * dz * (1 - dw) * table14 +
dx * dy * dz * dw * table15;
if (nw == 9) {
+ // interpolating metal species
dz = 0.0;
} else {
dw = 0.0;
@@ -475,7 +479,8 @@ __attribute__((always_inline)) INLINE float interpol_4d(
/*
* @brief Interpolates 2d EAGLE table over one of the dimensions,
- * producing 1d table
+ * producing 1d table. May be useful if need to use bisection
+ * method often with lots of iterations.
*
* @param p Particle structure
* @param cooling Cooling data structure
@@ -635,14 +640,6 @@ construct_1d_print_table_from_3d_elements(
result_table[i + array_size*j] += ((1 - d_x) * table[index[0]] +
d_x * table[index[1]])*abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * table[index[0]],
- (1 - d_x) * table[index[0]] +
- d_x * table[index[1]]);
-#endif
}
}
}
@@ -701,14 +698,6 @@ construct_1d_table_from_3d_elements(
result_table[i] += ((1 - d_x) * table[index[0]] +
d_x * table[index[1]])*abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * table[index[0]],
- (1 - d_x) * table[index[0]] +
- d_x * table[index[1]]);
-#endif
}
}
}
@@ -788,15 +777,6 @@ __attribute__((always_inline)) INLINE void construct_1d_table_from_4d(
d_x * (1 - d_y) * d_w * table[index[5]] +
d_x * d_y * (1 - d_w) * table[index[6]] +
d_x * d_y * d_w * table[index[7]];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 2 i dz dnh dHe table values %d %.5e %.5e %.5e %.5e "
- "%.5e %.5e %.5e %.5e %.5e %.5e %.5e \n",
- i, d_x, d_y, d_w, table[index[0]], table[index[1]],
- table[index[2]], table[index[3]], table[index[4]], table[index[5]],
- table[index[6]], table[index[7]]);
-#endif
}
}
@@ -861,21 +841,6 @@ construct_1d_print_table_from_4d_elements(
(1 - d_x) * d_y * table[index[1]] +
d_x * (1 - d_y) * table[index[2]] +
d_x * d_y * table[index[3]]) * abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * (1 - d_y) * table[index[0]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]);
-#endif
}
}
}
@@ -940,21 +905,6 @@ construct_1d_table_from_4d_elements(
(1 - d_x) * d_y * table[index[1]] +
d_x * (1 - d_y) * table[index[2]] +
d_x * d_y * table[index[3]]) * abundance_ratio[j+2];
-#if SWIFT_DEBUG_CHECKS
- if (isnan(result_table[i]))
- printf(
- "Eagle cooling.h 3 i partial sums %d %.5e %.5e %.5e %.5e %.5e \n",
- i, result_table[i], (1 - d_x) * (1 - d_y) * table[index[0]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]],
- (1 - d_x) * (1 - d_y) * table[index[0]] +
- (1 - d_x) * d_y * table[index[1]] +
- d_x * (1 - d_y) * table[index[2]] +
- d_x * d_y * table[index[3]]);
-#endif
}
}
}
@@ -996,20 +946,16 @@ eagle_convert_temp_to_u_1d_table(double temp, float *temperature_table,
* @param d_z Redshift offset
* @param d_n_h Hydrogen number density offset
* @param d_He Helium fraction offset
- * @param p Particle data structure
* @param cooling Cooling data structure
* @param cosmo Cosmology data structure
- * @param internal_const Physical constants 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 part *restrict p,
const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const) {
+ const struct cosmology *restrict cosmo) {
int u_i;
float d_u, logT;
@@ -1044,18 +990,12 @@ eagle_convert_u_to_temp(
* @param log_10_u Log base 10 of internal energy
* @param delta_u Pointer to size of internal energy cell
* @param temperature_table 1d array of temperatures
- * @param p Particle data structure
* @param cooling Cooling data structure
- * @param cosmo Cosmology data structure
- * @param internal_const Physical constants data structure
*/
__attribute__((always_inline)) INLINE double
eagle_convert_u_to_temp_1d_table(
double log_10_u, float *delta_u, double *temperature_table,
- const struct part *restrict p,
- const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
- const struct phys_const *internal_const) {
+ const struct cooling_function_data *restrict cooling) {
int u_i;
float d_u, logT;
@@ -1080,7 +1020,6 @@ eagle_convert_u_to_temp_1d_table(
* @param He_i Helium fraction index
* @param d_He Helium fraction offset
* @param phys_const Physical constants data structure
- * @param us Units data structure
* @param cosmo Cosmology data structure
* @param cooling Cooling data structure
* @param p Particle data structure
@@ -1110,7 +1049,6 @@ __attribute__((always_inline)) INLINE void construct_1d_tables(
int z_index, float dz, int n_h_i, float d_n_h,
int He_i, float d_He,
const struct phys_const *restrict phys_const,
- const struct unit_system *restrict us,
const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct part *restrict p,
diff --git a/src/cooling/EAGLE/interpolate.h b/src/cooling/EAGLE/interpolate.h
index 4d317bc668769b94aa255c08fbedc517a0a6230b..0392f234a601c9092bc32d6b8f6ee8e8ef735cd6 100644
--- a/src/cooling/EAGLE/interpolate.h
+++ b/src/cooling/EAGLE/interpolate.h
@@ -128,24 +128,18 @@ eagle_convert_u_to_temp(
double, float *,
int, int, int,
float, float, float,
- const struct part *restrict,
const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *);
+ const struct cosmology *restrict);
double
eagle_convert_u_to_temp_1d_table(
double, float *, double *,
- const struct part *restrict,
- const struct cooling_function_data *restrict,
- const struct cosmology *restrict,
- const struct phys_const *);
+ const struct cooling_function_data *restrict);
void construct_1d_tables(
int, float, int, float,
int, float,
const struct phys_const *restrict,
- const struct unit_system *restrict,
const struct cosmology *restrict,
const struct cooling_function_data *restrict,
const struct part *restrict,
diff --git a/src/engine.h b/src/engine.h
index 2a5a58793620e765755156d9cbcbed211312b2a1..18045d2e0c4df906e300304b4efa3833e160ec93 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -50,12 +50,6 @@
#include "task.h"
#include "units.h"
-// cooling counters
-int n_eagle_cooling_rate_calls_1;
-int n_eagle_cooling_rate_calls_2;
-int n_eagle_cooling_rate_calls_3;
-int n_eagle_cooling_rate_calls_4;
-
/**
* @brief The different policies the #engine can follow.
*/