Skip to content
Snippets Groups Projects
Commit a8e75a31 authored by Alexei Borissov's avatar Alexei Borissov Committed by Alexei Borissov
Browse files

update to automated cooling box test, general tidy

parent 4b9ae14f
No related branches found
No related tags found
1 merge request!593Eagle cooling
Hide whitespace changes
Inline Side-by-side
Showing
with 1077 additions and 2702 deletions
examples/CoolingBox/analytical_test.py
+ 45
− 66
View file @ a8e75a31
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]))
examples/CoolingBox/makeIC.py
+ 2
− 2
View file @ a8e75a31
......@@ -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"
......
examples/CoolingBox/run.sh
+ 1
− 1
View file @ a8e75a31
......@@ -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
examples/CoolingBox/test_script.sh
+ 17
− 15
View file @ a8e75a31
......@@ -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
......
examples/CoolingRates/plots.py 0 → 100644
+ 58
− 0
View file @ a8e75a31
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()
examples/CoolingRates/plots_nh.py 0 → 100644
+ 47
− 0
View file @ a8e75a31
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()
examples/CoolingRates/testCooling.c
+ 430
− 38
View file @ a8e75a31
......@@ -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;
......
examples/CoolingRates/testCooling.yml
+ 10
− 10
View file @ a8e75a31
......@@ -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
......
src/chemistry/EAGLE/chemistry.h
+ 0
− 2
View file @ a8e75a31
......@@ -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;
}
......
src/cooling/EAGLE/cooling.c
+ 66
− 135
View file @ a8e75a31
......@@ -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));
......
src/cooling/EAGLE/cooling.h
+ 118
− 2262
View file @ a8e75a31
......@@ -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 */
src/cooling/EAGLE/cooling_struct.h
+ 2
− 0
View file @ a8e75a31
......@@ -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
......
src/cooling/EAGLE/eagle_cool_tables.c
+ 249
− 65
View file @ a8e75a31
......@@ -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
src/cooling/EAGLE/interpolate.c
+ 30
− 92
View file @ a8e75a31
......@@ -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,
......
src/cooling/EAGLE/interpolate.h
+ 2
− 8
View file @ a8e75a31
......@@ -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,
......
src/engine.h
+ 0
− 6
View file @ a8e75a31
......@@ -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.
*/
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment