diff --git a/Makefile.am b/Makefile.am
index fb4eb5f6d6b63a7d0e034e0a3202ac61066e6e25..f92c9f27b046ec067cbc0d9b89861d4bf78b0f07 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -19,7 +19,7 @@
ACLOCAL_AMFLAGS = -I m4
# Show the way...
-SUBDIRS = src examples doc tests
+SUBDIRS = src examples doc tests examples/CoolingTest_Alexei
# Non-standard files that should be part of the distribution.
EXTRA_DIST = INSTALL.swift .clang-format format.sh
diff --git a/configure.ac b/configure.ac
index 24afb459d06cd71b855ac5598e10b8885c57bc36..bfab239722a1a97be78744f8ce6fb558facce5f4 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1556,7 +1556,7 @@ DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
# Handle .in files.
-AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
+AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/CoolingTest_Alexei/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([tests/testReading.sh], [chmod +x tests/testReading.sh])
AC_CONFIG_FILES([tests/testActivePair.sh], [chmod +x tests/testActivePair.sh])
AC_CONFIG_FILES([tests/test27cells.sh], [chmod +x tests/test27cells.sh])
diff --git a/examples/CoolingBox/analytical_test.py b/examples/CoolingBox/analytical_test.py
index 281de2b9b90f65e4c7e483778c197f71261b923f..62c33d0d87cf7af8beee92d6f6955befc3f53973 100644
--- a/examples/CoolingBox/analytical_test.py
+++ b/examples/CoolingBox/analytical_test.py
@@ -46,8 +46,10 @@ 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*rho/(0.59*m_p))
+ temp = pressure/(k_b*n)
+ #temp = pressure/(k_b*rho/(0.59*m_p))
return temp
# File containing the total energy
@@ -55,6 +57,7 @@ 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
@@ -65,7 +68,9 @@ 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
@@ -84,9 +89,10 @@ 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
-nsnap = 110
+nsnap = 501
npart = 4096
u_snapshots_cgs = zeros(nsnap)
u_part_snapshots_cgs = zeros((nsnap,npart))
@@ -96,6 +102,17 @@ for i in range(nsnap):
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
+
# calculate analytic solution. Since cooling rate is constant,
# temperature decrease in linear in time.
# read Lambda and temperatures from table
@@ -108,44 +125,49 @@ for line in file_in:
cooling_rate.append(-float(data[1]))
tfinal = t_snapshots_cgs[nsnap-1]
-nt = 1e4
+nt = 1e5
dt = tfinal/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_sol[0] = u
for j in range(int(nt-1)):
t_sol[j+1] = t_sol[j] + dt
- Lambda_net = interpol_lambda(temperature,cooling_rate,temp_sol[j])
- #u_next = (u*m_p - Lambda_net*rho/(0.59*m_p)*dt)/m_p
- nH = 0.7*rho/(m_p)
- #nH = 9.125e-5
+ 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
- #print(u_next, u, Lambda_net,rho/(0.59*m_p),dt)
- #print(nH**2/rho*Lambda_net, dt)
temp_sol[j+1] = convert_u_to_temp_sol(u_next,rho)
- lambda_sol[j] = Lambda_net
+ u_sol[j+1] = u_next
+ #lambda_sol[j] = Lambda_net
u = u_next
+
+print(u, Lambda_net)
mean_temp = np.zeros(nsnap)
+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)
-p = figure()
-p1, = plot(t_sol,temp_sol,linewidth = 0.5,color = 'k',label = 'analytical')
-p2, = plot(t_snapshots_cgs,mean_temp,linewidth = 0.5,color = 'r',label = 'swift mean')
+ 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)
+
+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')
l = legend(handles = [p1,p2])
-plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
-xlabel("${\\rm{Time~[s]}}$", labelpad=0)
-ylabel("${\\rm{Temperature~[K]}}$")
+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])
savefig("energy.png", dpi=200)
-#p = figure()
-#p1, = loglog(temp_sol,lambda_sol,linewidth = 0.5,color = 'k',label = 'analytical')
-#xlabel("${\\rm{Temperature~[K]}}$")
-#ylabel("${\\rm{Cooling~rate}}$")
-
-#savefig("cooling.png", dpi=200)
+print('Final internal energy ode, swift, error ' + str(u_sol[int(nt)-1]) + ' ' + str(mean_u[nsnap-1]) + ' ' + str( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])))
+#print(u_sol[int(nt)-1], mean_u[nsnap-1], (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1]))
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 4c376663bf1987203334d1199295e04c65c7c13d..aa0597c05fecabf007af0b4cda84bbf4f324c2b0 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -27,8 +27,8 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec
-rho = 3.2e6 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 4.5e7 # Pressure in code units (at 10^6K)
+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)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
diff --git a/examples/CoolingTest_Alexei/Makefile.am b/examples/CoolingTest_Alexei/Makefile.am
new file mode 100644
index 0000000000000000000000000000000000000000..e878a3792a84f7ace5f3b60fcc5875169145731b
--- /dev/null
+++ b/examples/CoolingTest_Alexei/Makefile.am
@@ -0,0 +1,32 @@
+# tHIS FIle is part of SWIFT.
+# Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+# Add the source directory and the non-standard paths to the included library headers to CFLAGS
+AM_CFLAGS = -I$(top_srcdir)/src $(HDF5_CPPFLAGS)
+
+AM_LDFLAGS = $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS) $(TCMALLOC_LIBS) $(JEMALLOC_LIBS) $(TBBMALLOC_LIBS) $(GRACKLE_LIBS) $(GSL_LIBS) $(PROFILER_LIBS)
+
+# Extra libraries.
+EXTRA_LIBS = $(HDF5_LIBS)
+
+# Programs.
+bin_PROGRAMS = testCooling
+
+# Sources
+testCooling_SOURCES = testCooling.c
+testCooling_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
+testCooling_LDADD = ../../src/.libs/libswiftsim.a $(EXTRA_LIBS)
+
diff --git a/examples/CoolingTest_Alexei/plots.py b/examples/CoolingTest_Alexei/plots.py
new file mode 100644
index 0000000000000000000000000000000000000000..969e6de713a4c421cc14d5f88ef68bafaab59a41
--- /dev/null
+++ b/examples/CoolingTest_Alexei/plots.py
@@ -0,0 +1,58 @@
+import matplotlib.pyplot as plt
+
+elements = 11
+temperature = []
+cooling_rate = [[] for i in range(elements+1)]
+u = []
+fu = []
+length = 0
+file_in = open('cooling_output.dat', 'r')
+for line in file_in:
+ data = line.split()
+ temperature.append(float(data[0]))
+ cooling_rate[0].append(-float(data[1]))
+
+file_in.close()
+
+for elem in range(elements):
+ file_in = open('cooling_element_'+str(elem)+'.dat','r')
+ for line in file_in:
+ data = line.split()
+ cooling_rate[elem+1].append(-float(data[0]))
+ file_in.close()
+
+#file_in = open('newton_output.dat', 'r')
+#for line in file_in:
+# data = line.split()
+# u.append(float(data[0]))
+# fu.append(float(data[1]))
+#
+#file_in.close()
+
+#p0, = plt.plot(u,fu, color = 'k')
+#p1, = plt.plot(u,[0 for x in u], color = 'r')
+#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
+#plt.xlim([1e13,2.0e14])
+#plt.ylim([0e13,2e14])
+#plt.xlabel('u')
+#plt.ylabel('f(u)')
+#plt.show()
+
+p0, = plt.loglog(temperature, cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(temperature, cooling_rate[1], linewidth = 0.5, color = 'k', linestyle = '--', label = 'H + He')
+p2, = plt.loglog(temperature, cooling_rate[3], linewidth = 0.5, color = 'b', label = 'Carbon')
+p3, = plt.loglog(temperature, cooling_rate[4], linewidth = 0.5, color = 'g', label = 'Nitrogen')
+p4, = plt.loglog(temperature, cooling_rate[5], linewidth = 0.5, color = 'r', label = 'Oxygen')
+p5, = plt.loglog(temperature, cooling_rate[6], linewidth = 0.5, color = 'c', label = 'Neon')
+p6, = plt.loglog(temperature, cooling_rate[7], linewidth = 0.5, color = 'm', label = 'Magnesium')
+p7, = plt.loglog(temperature, cooling_rate[8], linewidth = 0.5, color = 'y', label = 'Silicon')
+p8, = plt.loglog(temperature, cooling_rate[9], linewidth = 0.5, color = 'lightgray', label = 'Sulphur')
+p9, = plt.loglog(temperature, cooling_rate[10], linewidth = 0.5, color = 'olive', label = 'Calcium')
+p10, = plt.loglog(temperature, cooling_rate[11], linewidth = 0.5, color = 'saddlebrown', label = 'Iron')
+#plt.ylim([1e-24,1e-21])
+#plt.xlim([1e2,1e8])
+plt.xlabel('Temperature (K)')
+plt.ylabel('Cooling rate (eagle units...)')
+plt.legend(handles = [p0,p1,p2,p3,p4,p5,p6,p7,p8,p9,p10])
+#plt.legend(handles = [p0,p2,p3,p4,p5,p6,p7,p8,p9])
+plt.show()
diff --git a/examples/CoolingTest_Alexei/plots_nh.py b/examples/CoolingTest_Alexei/plots_nh.py
new file mode 100644
index 0000000000000000000000000000000000000000..ab3300750e2497ec8026dbe2d5c6c7d358f87831
--- /dev/null
+++ b/examples/CoolingTest_Alexei/plots_nh.py
@@ -0,0 +1,47 @@
+import matplotlib.pyplot as plt
+
+elements = 6
+temperature = [[] for i in range(elements+1)]
+cooling_rate = [[] for i in range(elements+1)]
+u = []
+fu = []
+length = 0
+
+for elem in range(elements):
+ file_in = open('cooling_output_'+str(elem)+'.dat','r')
+ for line in file_in:
+ data = line.split()
+ temperature[elem+1].append(float(data[0]))
+ cooling_rate[elem+1].append(-float(data[1]))
+ file_in.close()
+
+#file_in = open('newton_output.dat', 'r')
+#for line in file_in:
+# data = line.split()
+# u.append(float(data[0]))
+# fu.append(float(data[1]))
+#
+#file_in.close()
+
+#p0, = plt.plot(u,fu, color = 'k')
+#p1, = plt.plot(u,[0 for x in u], color = 'r')
+#p1, = plt.loglog(u,[-x for x in fu], color = 'r')
+#plt.xlim([1e13,2.0e14])
+#plt.ylim([0e13,2e14])
+#plt.xlabel('u')
+#plt.ylabel('f(u)')
+#plt.show()
+
+#p0, = plt.loglog(temperature[0], cooling_rate[0], linewidth = 0.5, color = 'k', label = 'Total')
+p1, = plt.loglog(temperature[1], cooling_rate[1], linewidth = 0.5, color = 'k', label = 'nh = 10^0')
+p2, = plt.loglog(temperature[2], cooling_rate[2], linewidth = 0.5, color = 'b', label = 'nh = 10^-1')
+p3, = plt.loglog(temperature[3], cooling_rate[3], linewidth = 0.5, color = 'g', label = 'nh = 10^-2')
+p4, = plt.loglog(temperature[4], cooling_rate[4], linewidth = 0.5, color = 'r', label = 'nh = 10^-3')
+p5, = plt.loglog(temperature[5], cooling_rate[5], linewidth = 0.5, color = 'c', label = 'nh = 10^-4')
+p6, = plt.loglog(temperature[6], cooling_rate[6], linewidth = 0.5, color = 'm', label = 'nh = 10^-5')
+plt.ylim([1e-24,1e-21])
+plt.xlim([1e4,1e8])
+plt.legend(handles = [p1,p2,p3,p4,p5,p6])
+plt.xlabel('Temperature (K)')
+plt.ylabel('Cooling rate (eagle units...)')
+plt.show()
diff --git a/examples/CoolingTest_Alexei/testCooling b/examples/CoolingTest_Alexei/testCooling
new file mode 100755
index 0000000000000000000000000000000000000000..20183de27bda2a0ecb577a27ab1c103024e98178
Binary files /dev/null and b/examples/CoolingTest_Alexei/testCooling differ
diff --git a/examples/CoolingTest_Alexei/testCooling.c b/examples/CoolingTest_Alexei/testCooling.c
new file mode 100644
index 0000000000000000000000000000000000000000..273f7a26812d7d898f010c3a2eeea38a289b362b
--- /dev/null
+++ b/examples/CoolingTest_Alexei/testCooling.c
@@ -0,0 +1,416 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "cooling.h"
+#include "hydro.h"
+#include "physical_constants.h"
+#include "swift.h"
+#include "units.h"
+
+enum table_index {
+ EAGLE_Hydrogen = 0,
+ EAGLE_Helium,
+ EAGLE_Carbon,
+ EAGLE_Nitrogen,
+ EAGLE_Oxygen,
+ EAGLE_Neon,
+ EAGLE_Magnesium,
+ EAGLE_Silicon,
+ EAGLE_Iron
+};
+
+void set_quantities(struct part *restrict p,
+ const struct unit_system *restrict us,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *restrict internal_const,
+ float nh,
+ double u){
+
+ const float gamma = 5.0/3.0;
+ double hydrogen_number_density = nh*pow(units_cgs_conversion_factor(us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo->z,3)));
+ p->rho = hydrogen_number_density*internal_const->const_proton_mass*
+ (1.0+p->chemistry_data.metal_mass_fraction[EAGLE_Helium]/p->chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+
+ float scale_factor = 1.0/(1.0+cosmo->z);
+ float pressure = (u*pow(scale_factor,3))/cooling->internal_energy_scale*p->rho*(gamma -1.0);
+ p->entropy = pressure/(pow(p->rho,gamma));
+}
+
+void construct_1d_tables_test(const struct part *restrict p,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo,
+ const struct phys_const *restrict internal_const,
+ double *temp_table,
+ double *H_plus_He_heat_table,
+ double *H_plus_He_electron_abundance_table,
+ double *element_electron_abundance_table,
+ double *element_print_cooling_table,
+ double *element_cooling_table,
+ float *abundance_ratio,
+ float *ub, float *lb){
+
+ float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
+ float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
+ float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale;
+
+ int z_index,He_i,n_h_i;
+ float dz,d_He,d_n_h;
+
+ get_redshift_index(cosmo->z, &z_index, &dz, cooling);
+ get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
+ get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
+ printf("testCooling 1d z_index dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n", z_index,dz, n_h_i, d_n_h, He_i, d_He);
+
+ construct_1d_table_from_4d(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.temperature,
+ z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
+ cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, temp_table, ub, lb);
+ construct_1d_table_from_4d(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.H_plus_He_heating, z_index, dz, cooling->N_Redshifts, n_h_i, d_n_h,
+ cooling->N_nH, He_i, d_He, cooling->N_He, cooling->N_Temp, H_plus_He_heat_table, ub, lb);
+ construct_1d_table_from_4d(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
+ dz, cooling->N_Redshifts, n_h_i, d_n_h, cooling->N_nH, He_i, d_He,
+ cooling->N_He, cooling->N_Temp, H_plus_He_electron_abundance_table, ub, lb);
+ construct_1d_table_from_4d_elements(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
+ n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_cooling_table,abundance_ratio, ub, lb);
+ construct_1d_print_table_from_4d_elements(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.metal_heating, z_index, dz, cooling->N_Redshifts,
+ n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_print_cooling_table,abundance_ratio, ub, lb);
+ construct_1d_table_from_3d(
+ p, cooling, cosmo, internal_const,
+ cooling->table.element_cooling.electron_abundance, z_index, dz, cooling->N_Redshifts,
+ n_h_i, d_n_h, cooling->N_nH, cooling->N_Temp, element_electron_abundance_table, ub, lb);
+
+}
+
+void compare_dlambda_du(
+ const struct unit_system *restrict us,
+ struct part *restrict p,
+ const struct xpart *restrict xp,
+ const struct phys_const *restrict internal_const,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo){
+
+ double H_plus_He_heat_table[176];
+ double H_plus_He_electron_abundance_table[176];
+ double temp_table[176];
+ double element_cooling_table[176];
+ double element_print_cooling_table[9 * 176];
+ double element_electron_abundance_table[176];
+ double rate_element_table[11];
+ float *abundance_ratio, cooling_du_dt1, cooling_du_dt2;
+ double dlambda_du1, dlambda_du2;
+
+ abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ abundance_ratio_to_solar(p, cooling, abundance_ratio);
+
+ float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
+ float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
+
+ int z_index,He_i,n_h_i;
+ float dz,d_He,d_n_h;
+ get_redshift_index(cosmo->z, &z_index, &dz, cooling);
+ get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
+
+ float nh = 1.0e-1, u = 1.0e9;
+ set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
+
+ const float log_10_e = 0.43429448190325182765;
+ float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
+ float lower_bound = cooling->Temp[0]/log_10_e;
+ construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
+ H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ element_electron_abundance_table, element_print_cooling_table,
+ element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
+
+ float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale;
+ get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
+ printf("testCooling.c 4d z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
+ int nt = 10;
+ for(int i = 0; i < nt; i++){
+ u = pow(10.0,9 + i);
+ set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
+ cooling_du_dt1 = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du1,H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_print_cooling_table,
+ element_electron_abundance_table,
+ temp_table,
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ p,cooling,cosmo,internal_const,rate_element_table,
+ abundance_ratio);
+ cooling_du_dt2 = eagle_metal_cooling_rate(log10(u),
+ &dlambda_du2,z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ p,cooling,cosmo,internal_const,NULL,
+ abundance_ratio);
+ printf("u du_dt_1d du_dt_4d dlambda_du_1d dlambda_du_4d %.5e %.5e %.5e %.5e %.5e\n",u,cooling_du_dt1,cooling_du_dt2,dlambda_du1,dlambda_du2);
+ }
+}
+
+void compare_temp(
+ const struct unit_system *restrict us,
+ struct part *restrict p,
+ const struct xpart *restrict xp,
+ const struct phys_const *restrict internal_const,
+ const struct cooling_function_data *restrict cooling,
+ const struct cosmology *restrict cosmo){
+
+ double H_plus_He_heat_table[176];
+ double H_plus_He_electron_abundance_table[176];
+ double temp_table[176];
+ double element_cooling_table[176];
+ double element_print_cooling_table[9 * 176];
+ double element_electron_abundance_table[176];
+ float *abundance_ratio;
+
+ abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ abundance_ratio_to_solar(p, cooling, abundance_ratio);
+
+ float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
+ float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
+
+ int z_index,He_i,n_h_i;
+ float dz,d_He,d_n_h;
+ get_redshift_index(cosmo->z, &z_index, &dz, cooling);
+ get_index_1d(cooling->HeFrac, cooling->N_He, HeFrac, &He_i, &d_He);
+
+ float nh = 1.0e-1, u = 1.0e9;
+ set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
+ const float log_10_e = 0.43429448190325182765;
+ float upper_bound = cooling->Temp[cooling->N_Temp-1]/log_10_e;
+ float lower_bound = cooling->Temp[0]/log_10_e;
+
+ construct_1d_tables_test(p, cooling, cosmo, internal_const, temp_table,
+ H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ element_electron_abundance_table, element_print_cooling_table,
+ element_cooling_table, abundance_ratio, &upper_bound, &lower_bound);
+ float T1d, T4d, delta_u;
+ float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
+ internal_const) *
+ cooling->number_density_scale;
+ get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
+ printf("testCooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e \n", z_index, dz, n_h_i, d_n_h, He_i, d_He);
+ int nt = 10;
+ for(int i = 0; i < nt; i++){
+ u = pow(10.0,9 + i);
+ set_quantities(p, us, cooling, cosmo, internal_const, nh, u);
+ T1d = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
+ p,cooling, cosmo, internal_const);
+ T4d = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
+ dz, d_n_h, d_He, p,cooling, cosmo, internal_const);
+ printf("u T1d T4d %.5e %.5e %.5e\n",u,pow(10.0,T1d),pow(10.0,T4d));
+ }
+}
+
+int main(int argc, char **argv) {
+ /* Declare relevant structs */
+ struct swift_params *params = malloc(sizeof(struct swift_params));
+ struct unit_system us;
+ struct chemistry_global_data chem_data;
+ struct part p;
+ struct xpart xp;
+ struct phys_const internal_const;
+ struct cooling_function_data cooling;
+ struct cosmology cosmo;
+ char *parametersFileName = "./testCooling.yml";
+
+ int tables = 0;
+
+ /* Read the parameter file */
+ if (params == NULL) error("Error allocating memory for the parameter file.");
+ message("Reading runtime parameters from file '%s'", parametersFileName);
+ parser_read_file(parametersFileName, params);
+
+ /* And dump the parameters as used. */
+ parser_write_params_to_file(params, "used_parameters.yml");
+
+ /* Init units */
+ units_init_from_params(&us, params, "InternalUnitSystem");
+ phys_const_init(&us, params, &internal_const);
+
+ /* Init chemistry */
+ chemistry_init(params, &us, &internal_const, &chem_data);
+ chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
+ chemistry_print(&chem_data);
+
+ /* Init cosmology */
+ cosmology_init(params, &us, &internal_const, &cosmo);
+ cosmology_print(&cosmo);
+
+ /* Init cooling */
+ cooling_init(params, &us, &internal_const, &cooling);
+ cooling_print(&cooling);
+
+ /* Calculate abundance ratios */
+ float *abundance_ratio;
+ abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
+
+ /* Declare 1D tables */
+ double H_plus_He_heat_table[176];
+ double H_plus_He_electron_abundance_table[176];
+ double temp_table[176];
+ double element_cooling_table[176];
+ double element_print_cooling_table[9 * 176];
+ double element_electron_abundance_table[176];
+
+ float nh;
+
+ float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
+ float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
+ int z_index,He_i,n_h_i;
+ float dz,d_He,d_n_h;
+ get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
+ get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
+
+ const float log_10_e = 0.43429448190325182765;
+ float upper_bound = cooling.Temp[cooling.N_Temp-1]/log_10_e;
+ float lower_bound = cooling.Temp[0]/log_10_e;
+
+ /* Loop over densities */
+ int nt = 250, n_nh = 6;
+ double u = pow(10.0,10);
+ if (argc == 1 || strcmp(argv[1], "nh") == 0){
+ for(int i = 0; i < n_nh; i++){
+ /* Open files */
+ char output_filename[21];
+ sprintf(output_filename, "%s%d%s", "cooling_output_", i, ".dat");
+ FILE *output_file = fopen(output_filename, "w");
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+
+ nh = pow(10.0,-i);
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
+ construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
+ temp_table, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_electron_abundance_table,
+ element_print_cooling_table,
+ element_cooling_table,
+ abundance_ratio, &upper_bound, &lower_bound);
+ get_index_1d(cooling.nH, cooling.N_nH, log10(nh), &n_h_i, &d_n_h);
+
+ for(int j = 0; j < nt; j++){
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,11.0 + j*8.0/nt));
+ u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
+ double dlambda_du;
+ float delta_u, cooling_du_dt, logT;
+ if (tables == 1) {
+ cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u),&dlambda_du,H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_print_cooling_table,
+ element_electron_abundance_table,
+ temp_table,
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,NULL,
+ abundance_ratio);
+ logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
+ &p,&cooling, &cosmo, &internal_const);
+ } else {
+ cooling_du_dt = eagle_metal_cooling_rate(log10(u),
+ &dlambda_du,z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,NULL,
+ abundance_ratio);
+ logT = eagle_convert_u_to_temp(log10(u), &delta_u, z_index, n_h_i, He_i,
+ dz, d_n_h, d_He, &p,&cooling, &cosmo, &internal_const);
+ }
+ float temperature_swift = pow(10.0,logT);
+
+ fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
+ }
+ fclose(output_file);
+ }
+ }
+ if (argc == 1 || strcmp(argv[1],"metals") == 0) {
+ printf("here \n");
+ char output_filename[21];
+ sprintf(output_filename, "%s", "cooling_output.dat");
+ FILE *output_file = fopen(output_filename, "w");
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ nh = pow(10.0,0);
+ u = pow(10.0,14.0);
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, u);
+ float nh_swift = chemistry_get_number_density(&p, &cosmo, chemistry_element_H, &internal_const)*cooling.number_density_scale;
+ printf("hydrogen number density %.5e\n", nh_swift);
+ construct_1d_tables_test(&p, &cooling, &cosmo, &internal_const,
+ temp_table, H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_electron_abundance_table,
+ element_print_cooling_table,
+ element_cooling_table,
+ abundance_ratio, &upper_bound, &lower_bound);
+ float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+ &internal_const) *
+ cooling.number_density_scale;
+ get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
+ printf("testcooling.c z_i dz nh_i d_nh He_i d_He %d %.5e %d %.5e %d %.5e\n",z_index, dz, n_h_i, d_n_h, He_i, d_He);
+ for(int j = 0; j < nt; j++){
+ set_quantities(&p, &us, &cooling, &cosmo, &internal_const, nh, pow(10.0,10.0 + j*8.0/nt));
+ u = hydro_get_physical_internal_energy(&p,&cosmo)*cooling.internal_energy_scale;
+ float delta_u, cooling_du_dt, logT;
+ if (tables == 1) {
+ cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_print_cooling_table,
+ element_electron_abundance_table,
+ temp_table,
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,
+ abundance_ratio);
+ logT = eagle_convert_u_to_temp_1d_table(log10(u), &delta_u, temp_table,
+ &p,&cooling, &cosmo, &internal_const);
+ } else {
+ cooling_du_dt = eagle_print_metal_cooling_rate(
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,
+ abundance_ratio);
+ }
+ fprintf(output_file,"%.5e %.5e\n", u,cooling_du_dt);
+ }
+ fclose(output_file);
+ }
+
+ if (argc == 1 || strcmp(argv[1],"compare_temp") == 0) compare_temp(&us, &p, &xp, &internal_const, &cooling, &cosmo);
+ if (argc == 1 || strcmp(argv[1],"compare_dlambda_du") == 0) compare_dlambda_du(&us, &p, &xp, &internal_const, &cooling, &cosmo);
+
+ free(params);
+ return 0;
+}
+
+
diff --git a/examples/CoolingTest_Alexei/testCooling.yml b/examples/CoolingTest_Alexei/testCooling.yml
new file mode 100644
index 0000000000000000000000000000000000000000..8426042ac7f2e7d18f9fc98b978bfcdad429be4a
--- /dev/null
+++ b/examples/CoolingTest_Alexei/testCooling.yml
@@ -0,0 +1,94 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.989e43 # 10^10 M_sun in grams
+ UnitLength_in_cgs: 3.085678e24 # Mpc in centimeters
+ UnitVelocity_in_cgs: 1e5 # km/s in centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Cosmological parameters
+Cosmology:
+ h: 0.6777 # Reduced Hubble constant
+ #a_begin: 0.9090909 # Initial scale-factor of the simulation
+ a_begin: 1.0 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
+ Omega_m: 0.307 # Matter density parameter
+ Omega_lambda: 0.693 # Dark-energy density parameter
+ Omega_b: 0.0455 # Baryon density parameter
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+Scheduler:
+ max_top_level_cells: 15
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: eagle # Common part of the name of output files
+ scale_factor_first: 0.92 # Scale-factor of the first snaphot (cosmological run)
+ time_first: 0.01 # Time of the first output (non-cosmological run) (in internal units)
+ delta_time: 1.10 # Time difference between consecutive outputs (in internal units)
+ compression: 1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ scale_factor_first: 0.92 # Scale-factor of the first stat dump (cosmological run)
+ time_first: 0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+ delta_time: 1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./EAGLE_ICs_12.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
+
+EAGLEChemistry:
+ InitMetallicity: 0.014
+ InitAbundance_Hydrogen: 0.70649785
+ InitAbundance_Helium: 0.28055534
+ InitAbundance_Carbon: 2.0665436e-3
+ InitAbundance_Nitrogen: 8.3562563e-4
+ InitAbundance_Oxygen: 5.4926244e-3
+ InitAbundance_Neon: 1.4144605e-3
+ InitAbundance_Magnesium: 5.907064e-4
+ InitAbundance_Silicon: 6.825874e-4
+ InitAbundance_Iron: 1.1032152e-3
+ CalciumOverSilicon: 0.0941736
+ SulphurOverSilicon: 0.6054160
+ #InitMetallicity: 0.
+ #InitAbundance_Hydrogen: 0.752
+ #InitAbundance_Helium: 0.248
+ #InitAbundance_Carbon: 0.000
+ #InitAbundance_Nitrogen: 0.000
+ #InitAbundance_Oxygen: 0.000
+ #InitAbundance_Neon: 0.000
+ #InitAbundance_Magnesium: 0.000
+ #InitAbundance_Silicon: 0.000
+ #InitAbundance_Iron: 0.000
+ #CalciumOverSilicon: 0.0941736
+ #SulphurOverSilicon: 0.6054160
+
+EagleCooling:
+ filename: /cosma5/data/Eagle/BG_Tables/CoolingTables/
+ reionisation_redshift: 8.898
+ he_reion: 1
+ he_reion_z_center: 3.5
+ he_reion_z_sigma: 0.5
+ he_reion_ev_pH: 2.0
+
diff --git a/examples/CoolingTest_Alexei/testCooling_messy.c b/examples/CoolingTest_Alexei/testCooling_messy.c
new file mode 100644
index 0000000000000000000000000000000000000000..b5c00e5625c43e57e46c5f8426640b9f511bc58a
--- /dev/null
+++ b/examples/CoolingTest_Alexei/testCooling_messy.c
@@ -0,0 +1,354 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "cooling.h"
+#include "hydro.h"
+#include "physical_constants.h"
+#include "swift.h"
+#include "units.h"
+
+int main(int argc, char **argv) {
+
+ /* Read the parameter file */
+ struct swift_params *params = malloc(sizeof(struct swift_params));
+ struct unit_system us;
+ struct chemistry_global_data chem_data;
+ struct part p;
+ struct xpart xp;
+ struct phys_const internal_const;
+ struct cooling_function_data cooling;
+ struct cosmology cosmo;
+ char *parametersFileName = "./testCooling.yml";
+ enum table_index {
+ EAGLE_Hydrogen = 0,
+ EAGLE_Helium,
+ EAGLE_Carbon,
+ EAGLE_Nitrogen,
+ EAGLE_Oxygen,
+ EAGLE_Neon,
+ EAGLE_Magnesium,
+ EAGLE_Silicon,
+ EAGLE_Iron
+ };
+
+ if (params == NULL) error("Error allocating memory for the parameter file.");
+ message("Reading runtime parameters from file '%s'", parametersFileName);
+ parser_read_file(parametersFileName, params);
+
+ /* And dump the parameters as used. */
+ parser_write_params_to_file(params, "used_parameters.yml");
+
+ double UnitMass_in_cgs = 1.989e43;
+ double UnitLength_in_cgs = 3.085678e24;
+ double UnitVelocity_in_cgs = 1e5;
+ double UnitCurrent_in_cgs = 1;
+ double UnitTemp_in_cgs = 1;
+ units_init(&us, UnitMass_in_cgs, UnitLength_in_cgs, UnitVelocity_in_cgs, UnitCurrent_in_cgs, UnitTemp_in_cgs);
+ phys_const_init(&us, params, &internal_const);
+
+ double hydrogen_number_density_cgs = 1e-4;
+ double u, u_cgs, hydrogen_number_density, pressure, gamma, cooling_du_dt,
+ temperature_cgs, newton_func, u_ini_cgs, ratefact, dt_cgs;
+ u_cgs = 0;
+ cooling_du_dt = 0;
+ temperature_cgs = 0;
+ newton_func = 0;
+ // int n_t_i = 2000;
+ dt_cgs = 4.0e-8 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
+
+ char *output_filename = "cooling_output.dat";
+ FILE *output_file = fopen(output_filename, "w");
+ if (output_file == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ char *output_filename2 = "temperature_output.dat";
+ FILE *output_file2 = fopen(output_filename2, "w");
+ if (output_file2 == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ char *output_filename3 = "newton_output.dat";
+ FILE *output_file3 = fopen(output_filename3, "w");
+ if (output_file3 == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+
+ gamma = 5.0 / 3.0;
+
+ chemistry_init(params, &us, &internal_const, &chem_data);
+ chemistry_first_init_part(&internal_const, &us, &cosmo, &chem_data, &p, &xp);
+ chemistry_print(&chem_data);
+
+ cosmology_init(params, &us, &internal_const, &cosmo);
+ cosmology_print(&cosmo);
+
+ u = 1.0 * pow(10.0, 11) /
+ (units_cgs_conversion_factor(&us, UNIT_CONV_ENERGY) /
+ units_cgs_conversion_factor(&us, UNIT_CONV_MASS));
+ pressure = u * p.rho * (gamma - 1.0);
+ hydrogen_number_density =
+ hydrogen_number_density_cgs *
+ pow(units_cgs_conversion_factor(&us, UNIT_CONV_LENGTH), 3);
+ p.rho = hydrogen_number_density * internal_const.const_proton_mass *
+ (1.0 +
+ p.chemistry_data.metal_mass_fraction[EAGLE_Helium] /
+ p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ p.entropy = pressure / (pow(p.rho, gamma));
+
+ cooling_init(params, &us, &internal_const, &cooling);
+ cooling_print(&cooling);
+
+ float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
+ float HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
+ float inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+ &internal_const) *
+ cooling.number_density_scale;
+ ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+
+ float *abundance_ratio;
+ abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
+
+ // construct 1d table of cooling rates wrt temperature
+ double H_plus_He_heat_table[176];
+ double H_plus_He_electron_abundance_table[176];
+ double temp_table[176];
+ double element_cooling_table[176];
+ double element_print_cooling_table[9 * 176];
+ double element_electron_abundance_table[176];
+
+ int z_index,He_i,n_h_i;
+ float dz,d_He,d_n_h;
+ get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
+ get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
+ get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
+
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.temperature,
+ z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
+ cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
+ cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
+ dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
+ cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
+ construct_1d_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
+ construct_1d_print_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
+ construct_1d_table_from_3d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
+
+
+ //
+ // printf("cooling table values \n");
+ // for( int j=0; j < 176; j++) {
+ // printf(" %.5e", H_plus_He_heat_table[j]+element_cooling_table[j]
+ u_ini_cgs = pow(10.0, 14);
+ float delta_u;
+
+
+ //Compute contributions to cooling rate from different metals
+ // int n_t_i = 100;
+ // for(int t_i = 0; t_i < n_t_i; t_i++){
+ // u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
+ // //u = u_cgs/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
+ // u = u_cgs/cooling.internal_energy_scale;
+ // hydrogen_number_density =
+ // hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
+ // p.rho =
+ // hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+ // pressure = u*p.rho*(gamma -1.0);
+ // p.entropy = pressure/(pow(p.rho,gamma));
+ // double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
+ // cooling.internal_energy_scale;
+
+ // cooling_du_dt = eagle_print_metal_cooling_rate_1d_table(H_plus_He_heat_table,
+ // H_plus_He_electron_abundance_table,
+ // element_print_cooling_table,
+ // element_electron_abundance_table,
+ // temp_table,
+ // z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ // &p,&cooling,&cosmo,&internal_const,
+ // abundance_ratio);
+ // float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
+ // &p,&cooling, &cosmo, &internal_const);
+ // float temperature_swift = pow(10.0,logT);
+
+ // fprintf(output_file,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
+ // fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
+ // temperature_cgs);
+
+ // newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
+ // fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
+ // //printf("u,u_ini,cooling_du_dt,ratefact,dt %.5e %.5e %.5e %.5e %.5e
+ // // \n",u_cgs,u_ini_cgs,cooling_du_dt,ratefact,dt_cgs);
+ //}
+
+ int n_t_i = 100;
+ for(int nh_i = 1; nh_i < 7; nh_i++){
+ char output_filename4[21];
+ sprintf(output_filename4, "%s%d%s", "cooling_output_", nh_i, ".dat");
+ FILE *output_file4 = fopen(output_filename4, "w");
+ if (output_file4 == NULL) {
+ printf("Error opening file!\n");
+ exit(1);
+ }
+ hydrogen_number_density = pow(10.0,-nh_i)*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3)*(1.0/(pow(1.0+cosmo.z,3)));
+ p.rho =
+ hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
+
+ XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
+ HeFrac = p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
+ inn_h = chemistry_get_number_density(&p, &cosmo, chemistry_element_H,
+ &internal_const) *
+ cooling.number_density_scale;
+ ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+ printf("test cooling hydrogen num dens inn_h %.5e %.5e\n", hydrogen_number_density*cooling.number_density_scale, inn_h);
+
+ abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
+
+ get_redshift_index(cosmo.z, &z_index, &dz, &cooling);
+ get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
+ get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
+
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.temperature,
+ z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
+ cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, temp_table);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_heating, z_index, dz, cooling.N_Redshifts, n_h_i, d_n_h,
+ cooling.N_nH, He_i, d_He, cooling.N_He, cooling.N_Temp, H_plus_He_heat_table);
+ construct_1d_table_from_4d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.H_plus_He_electron_abundance, z_index,
+ dz, cooling.N_Redshifts, n_h_i, d_n_h, cooling.N_nH, He_i, d_He,
+ cooling.N_He, cooling.N_Temp, H_plus_He_electron_abundance_table);
+ construct_1d_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_cooling_table,abundance_ratio);
+ construct_1d_print_table_from_4d_elements(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.metal_heating, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_print_cooling_table,abundance_ratio);
+ construct_1d_table_from_3d(
+ &p, &cooling, &cosmo, &internal_const,
+ cooling.table.element_cooling.electron_abundance, z_index, dz, cooling.N_Redshifts,
+ n_h_i, d_n_h, cooling.N_nH, cooling.N_Temp, element_electron_abundance_table);
+ for(int t_i = 0; t_i < n_t_i; t_i++){
+ u_cgs = pow(10.0,11 + t_i*6.0/n_t_i);
+ u = u_cgs/cooling.internal_energy_scale;
+ pressure = u*p.rho*(gamma -1.0);
+ p.entropy = pressure/(pow(p.rho,gamma));
+ double u_swift = hydro_get_physical_internal_energy(&p, &cosmo) *
+ cooling.internal_energy_scale;
+
+ double dlambda_du;
+ cooling_du_dt = eagle_metal_cooling_rate_1d_table(log10(u_swift),&dlambda_du,H_plus_He_heat_table,
+ H_plus_He_electron_abundance_table,
+ element_print_cooling_table,
+ element_electron_abundance_table,
+ temp_table,
+ z_index, dz, n_h_i, d_n_h, He_i, d_He,
+ &p,&cooling,&cosmo,&internal_const,NULL,
+ abundance_ratio);
+ float logT = eagle_convert_u_to_temp_1d_table(log10(u_swift), &delta_u, temp_table,
+ &p,&cooling, &cosmo, &internal_const);
+ float temperature_swift = pow(10.0,logT);
+
+ fprintf(output_file4,"%.5e %.5e\n", temperature_swift,cooling_du_dt);
+ fprintf(output_file2,"%.5e %.5e\n",u_swift*(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS)),
+ temperature_cgs);
+
+ newton_func = u_cgs - u_ini_cgs - cooling_du_dt*ratefact*dt_cgs;
+ fprintf(output_file3,"%.5e %.5e\n",u_cgs,newton_func);
+ }
+ fclose(output_file4);
+ }
+ fclose(output_file);
+ fclose(output_file2);
+ fclose(output_file3);
+
+ //double dLambdaNet_du, LambdaNet; //, LambdaNext;
+ //float x_init, u_eq = 2.0e12;
+ //for (int j = 5; j < 6; j++) {
+ // float u_ini = eagle_convert_temp_to_u_1d_table(pow(10.0, 0.5 * (j + 5)),
+ // temp_table, &p, &cooling,
+ // &cosmo, &internal_const),
+ // x, du;
+ // float dt = 2.0e-4 * units_cgs_conversion_factor(&us, UNIT_CONV_TIME);
+ // LambdaNet = eagle_cooling_rate_1d_table(
+ // u_ini, &dLambdaNet_du, H_plus_He_heat_table,
+ // H_plus_He_electron_abundance_table, element_cooling_table,
+ // element_electron_abundance_table, temp_table, &p, &cooling, &cosmo,
+ // &internal_const);
+ // float u_temp = u_ini + LambdaNet * ratefact * dt;
+ // /* RGB removed this **
+ // if (u_temp > 0) LambdaNext = eagle_cooling_rate_1d_table(u_temp,
+ // &dLambdaNet_du, H_plus_He_heat_table, H_plus_He_electron_abundance_table,
+ // element_cooling_table, element_electron_abundance_table, temp_table, &p,
+ // &cooling, &cosmo, &internal_const);
+ // if (fabs(LambdaNet - LambdaNext)/LambdaNet < 0.5) {
+ // u_temp = u_ini;
+ // } else {
+ // u_temp =
+ // eagle_convert_temp_to_u_1d_table(1.0e4,temp_table,&p,&cooling,&cosmo,&internal_const);
+ // } */
+ // if (u_temp > u_eq) {
+ // x_init = log(u_temp);
+ // } else {
+ // x_init = log(u_eq);
+ // }
+ // x = newton_iter(x_init, u_ini, H_plus_He_heat_table,
+ // H_plus_He_electron_abundance_table, element_cooling_table,
+ // element_electron_abundance_table, temp_table, &p, &cosmo,
+ // &cooling, &internal_const, dt);
+ // printf(
+ // "testing newton integration, u_ini, u %.5e %.5e, temperature initial, "
+ // "final %.5e %.5e\n",
+ // u_ini, exp(x),
+ // eagle_convert_u_to_temp_1d_table(u_ini, &du, temp_table, &p, &cooling,
+ // &cosmo, &internal_const),
+ // eagle_convert_u_to_temp_1d_table(exp(x), &du, temp_table, &p, &cooling,
+ // &cosmo, &internal_const));
+ //}
+
+ free(params);
+
+ return 0;
+}
diff --git a/src/cooling/EAGLE/.cooling.h.swp b/src/cooling/EAGLE/.cooling.h.swp
new file mode 100644
index 0000000000000000000000000000000000000000..dbf1ffd09a09fea79583fc056bc56c16c449684f
Binary files /dev/null and b/src/cooling/EAGLE/.cooling.h.swp differ
diff --git a/src/cooling/EAGLE/.eagle_cool_tables.h.swp b/src/cooling/EAGLE/.eagle_cool_tables.h.swp
new file mode 100644
index 0000000000000000000000000000000000000000..8e230bda82187c4455ad6dcefc873ee472ddd69a
Binary files /dev/null and b/src/cooling/EAGLE/.eagle_cool_tables.h.swp differ
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index b0530e517b6f68b0ac2c21e28ba6cd813c09a0ca..998bfdbc3230700ee4055ebe430f5a36be7daabb 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -528,6 +528,54 @@ __attribute__((always_inline)) INLINE float interpol_4d(
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
@@ -548,24 +596,37 @@ __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 z_i, float d_z,
- int n_h_i, float d_n_h, double *result_table) {
- int index[4];
+ 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--;
- for (int i = 0; i < cooling->N_Temp; i++) {
- index[0] = row_major_index_3d(z_i, n_h_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
- index[1] = row_major_index_3d(z_i, n_h_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
- index[2] = row_major_index_3d(z_i + 1, n_h_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
- index[3] = row_major_index_3d(z_i + 1, n_h_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_Temp);
-
- result_table[i] = (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]] +
- d_z * d_n_h * table[index[3]];
+ 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]];
}
}
@@ -598,74 +659,67 @@ 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 z_i, float d_z,
- int He_i, float d_He, int n_h_i, float d_n_h, double *result_table,
- double *x_max_grad, double u_ini, float dt) {
- int index[8];
- float x0, x1, y0, y1, old_grad, new_grad;
- float inn_h = chemistry_get_number_density(p, cosmo, chemistry_element_H,
- internal_const) *
- cooling->number_density_scale;
- float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
- float ratefact = inn_h * (XH / eagle_proton_mass_cgs);
-
- old_grad = 0.0;
+ 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--;
- for (int i = 0; i < cooling->N_Temp; i++) {
+ int index[8];
+ for (int i = index_low; i < index_high; i++) {
index[0] =
- row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i + 1, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i + 1, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i + 1, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[7] = row_major_index_4d(z_i + 1, n_h_i + 1, He_i + 1, i,
- cooling->N_Redshifts, cooling->N_nH,
- cooling->N_He, cooling->N_Temp);
-
- result_table[i] = (1 - d_z) * (1 - d_n_h) * (1 - d_He) * table[index[0]] +
- (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
- (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
- (1 - d_z) * d_n_h * d_He * table[index[3]] +
- d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
- d_z * (1 - d_n_h) * d_He * table[index[5]] +
- d_z * d_n_h * (1 - d_He) * table[index[6]] +
- d_z * d_n_h * d_He * table[index[7]];
+ 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_z, d_n_h, d_He, table[index[0]], table[index[1]],
+ 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
- if (i > 0 && dt > 0 && x_max_grad != NULL) {
- x0 = pow(10.0, cooling->Therm[i - 1]);
- x1 = pow(10.0, cooling->Therm[i]);
- y0 = x0 - u_ini - result_table[i - 1] * ratefact * dt;
- y1 = x1 - u_ini - result_table[i] * ratefact * dt;
- new_grad = log10(y1 / y0) / log10(x1 / x0);
- if (new_grad > old_grad) {
- *x_max_grad = 0.5 * (x0 + x1);
- old_grad = new_grad;
- }
- //printf("Eagle cooling.h id grad max_grad x_max_grad %llu %.5e %.5e %.5e\n", p->id, new_grad, old_grad, *x_max_grad);
- }
}
}
@@ -691,49 +745,63 @@ __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 z_i, float d_z,
- int He_i, float d_He, int n_h_i, float d_n_h, double *result_table) {
+ 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 = 0; i < cooling->N_Temp; i++) {
+ for (int i = index_low; i < index_high; i++) {
index[0] =
- row_major_index_4d(z_i, n_h_i, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i + 1, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i, n_h_i + 1, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i, He_i + 1, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ 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(z_i + 1, n_h_i + 1, He_i, i, cooling->N_Redshifts,
- cooling->N_nH, cooling->N_He, cooling->N_Temp);
- index[7] = row_major_index_4d(z_i + 1, n_h_i + 1, He_i + 1, i,
- cooling->N_Redshifts, cooling->N_nH,
- cooling->N_He, cooling->N_Temp);
-
- result_table[i] = (1 - d_z) * (1 - d_n_h) * (1 - d_He) * table[index[0]] +
- (1 - d_z) * (1 - d_n_h) * d_He * table[index[1]] +
- (1 - d_z) * d_n_h * (1 - d_He) * table[index[2]] +
- (1 - d_z) * d_n_h * d_He * table[index[3]] +
- d_z * (1 - d_n_h) * (1 - d_He) * table[index[4]] +
- d_z * (1 - d_n_h) * d_He * table[index[5]] +
- d_z * d_n_h * (1 - d_He) * table[index[6]] +
- d_z * d_n_h * d_He * table[index[7]];
+ 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_z, d_n_h, d_He, table[index[0]], table[index[1]],
+ 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
@@ -773,6 +841,88 @@ __attribute__((always_inline)) INLINE double eagle_helium_reionization_extraheat
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, j, y_i, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[1] = row_major_index_4d(x_i, j, y_i + 1, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[2] = row_major_index_4d(x_i + 1, j, y_i, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[3] = row_major_index_4d(x_i + 1, j, y_i + 1, i,
+ n_x, cooling->N_Elements,
+ n_y, array_size);
+
+ 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.
@@ -798,44 +948,123 @@ 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 z_i, float d_z,
- int n_h_i, float d_n_h, double *result_table) {
+ 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 = 0; i < cooling->N_Temp; i++) {
+ for (int i = index_low; i < index_high; i++) {
if (j == 0) result_table[i] = 0.0;
- index[0] = row_major_index_4d(z_i, j, n_h_i, i, cooling->N_Redshifts,
- cooling->N_Elements, cooling->N_nH,
- cooling->N_Temp);
- index[1] = row_major_index_4d(z_i, j, n_h_i + 1, i, cooling->N_Redshifts,
- cooling->N_Elements, cooling->N_nH,
- cooling->N_Temp);
- index[2] = row_major_index_4d(z_i + 1, j, n_h_i, i, cooling->N_Redshifts,
- cooling->N_Elements, cooling->N_nH,
- cooling->N_Temp);
- index[3] = row_major_index_4d(z_i + 1, j, n_h_i + 1, i,
- cooling->N_Redshifts, cooling->N_Elements,
- cooling->N_nH, cooling->N_Temp);
-
- result_table[i] += (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]] +
- d_z * d_n_h * table[index[3]];
+ index[0] = row_major_index_4d(x_i, j, y_i, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[1] = row_major_index_4d(x_i, j, y_i + 1, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[2] = row_major_index_4d(x_i + 1, j, y_i, i, n_x,
+ cooling->N_Elements, n_y,
+ array_size);
+ index[3] = row_major_index_4d(x_i + 1, j, y_i + 1, i,
+ n_x, cooling->N_Elements,
+ 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]]) * 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_z) * (1 - d_n_h) * table[index[0]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]],
- (1 - d_z) * (1 - d_n_h) * table[index[0]] +
- (1 - d_z) * d_n_h * table[index[1]] +
- d_z * (1 - d_n_h) * table[index[2]] +
- d_z * d_n_h * table[index[3]]);
+ 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
}
}
@@ -1003,6 +1232,41 @@ eagle_convert_temp_to_u_1d_table(double temp, float *temperature_table,
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;
+
+ 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);
+
+ *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
@@ -1015,7 +1279,7 @@ eagle_convert_temp_to_u_1d_table(double temp, float *temperature_table,
*/
__attribute__((always_inline)) INLINE static double
eagle_convert_u_to_temp_1d_table(
- double logu, float *delta_u, double *temperature_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,
@@ -1024,13 +1288,9 @@ eagle_convert_u_to_temp_1d_table(
int u_i;
float d_u, logT;//, T;
- get_index_1d(cooling->Therm, cooling->N_Temp, logu, &u_i, &d_u);
+ 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);
- //T = pow(10.0, logT);
-
- //if (u_i == 0 && d_u == 0) T *= u / pow(10.0, cooling->Therm[0]);
- if (u_i == 0 && d_u == 0) logT += logu - cooling->Therm[0];
*delta_u =
pow(10.0, cooling->Therm[u_i + 1]) - pow(10.0, cooling->Therm[u_i]);
@@ -1057,29 +1317,34 @@ eagle_convert_u_to_temp_1d_table(
* from each of the metals.
*/
__attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(
- double logu, double *temp_table, int z_index, float dz, int n_h_i, float d_n_h,
+ 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 *ratio_solar) {
- double T_gam, solar_electron_abundance;
+ 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;
+ double cooling_rate = 0.0, temp_lambda, temp_lambda1, temp_lambda2;
float du;
float h_plus_he_electron_abundance;
+ float h_plus_he_electron_abundance1;
+ float h_plus_he_electron_abundance2;
int i;
double temp;
int temp_i;
float d_temp;
- temp = eagle_convert_u_to_temp_1d_table(logu, &du, temp_table, p, cooling, cosmo,
- internal_const);
+ *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);
@@ -1113,11 +1378,29 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(
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_lambda1 = 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,
+ 0.0, cooling->N_Redshifts, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp);
+ temp_lambda2 = interpol_4d(cooling->table.element_cooling.H_plus_He_heating,
+ z_index, n_h_i, He_i, temp_i+1, dz, d_n_h, d_He,
+ 0.0, cooling->N_Redshifts, cooling->N_nH,
+ cooling->N_He, cooling->N_Temp);
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_abundance1 = 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, 0.0, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
+ h_plus_he_electron_abundance2 = interpol_4d(
+ cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
+ n_h_i, He_i, temp_i+1, dz, d_n_h, d_He, 0.0, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_He, cooling->N_Temp);
cooling_rate += temp_lambda;
+ //printf("Eagle cooling.h 4d H_plus_He_cooling electron abundance %.5e %.5e\n", temp_lambda, h_plus_he_electron_abundance);
+ *dlambda_du += (temp_lambda2 - temp_lambda1)/du;
if (element_lambda != NULL) element_lambda[0] = temp_lambda;
/* ------------------ */
@@ -1181,18 +1464,37 @@ __attribute__((always_inline)) INLINE static double eagle_metal_cooling_rate(
/* redshift tables */
solar_electron_abundance =
interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,
- n_h_i, dz, temp_i, d_n_h, d_temp, cooling->N_Redshifts,
+ n_h_i, temp_i, dz, d_n_h, d_temp, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ solar_electron_abundance1 =
+ interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,
+ n_h_i, temp_i, dz, d_n_h, 0.0, cooling->N_Redshifts,
+ cooling->N_nH, cooling->N_Temp);
+ solar_electron_abundance2 =
+ interpol_3d(cooling->table.element_cooling.electron_abundance, z_index,
+ n_h_i, temp_i+1, dz, d_n_h, 0.0, cooling->N_Redshifts,
cooling->N_nH, cooling->N_Temp);
for (i = 0; i < cooling->N_Elements; i++) {
- // WARNING: before doing proper runs need to
- // multiply by ratio of particle abundances to solar
temp_lambda =
interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i,
n_h_i, temp_i, dz, 0.0, d_n_h, d_temp, cooling->N_Redshifts,
cooling->N_Elements, cooling->N_nH, cooling->N_Temp) *
- (h_plus_he_electron_abundance / solar_electron_abundance);
+ (h_plus_he_electron_abundance / solar_electron_abundance)*
+ solar_ratio[i+2];
+ elem_cool1 = interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i,
+ n_h_i, temp_i, dz, 0.0, d_n_h, 0.0, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
+ elem_cool2 = interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i,
+ n_h_i, temp_i+1, dz, 0.0, d_n_h, 0.0, cooling->N_Redshifts,
+ cooling->N_Elements, cooling->N_nH, cooling->N_Temp);
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];
+ //printf("Eagle cooling.h 4d element cooling electron abundance, solar_ratio %d %.5e %.5e %.5e\n", i, temp_lambda, solar_electron_abundance, solar_ratio[i+2]);
if (element_lambda != NULL) element_lambda[i + 2] = temp_lambda;
}
@@ -1256,6 +1558,7 @@ eagle_metal_cooling_rate_1d_table(
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 */
@@ -1290,6 +1593,7 @@ eagle_metal_cooling_rate_1d_table(
*dlambda_du += (((double)H_plus_He_heat_table[temp_i + 1]) -
((double)H_plus_He_heat_table[temp_i])) /
((double)du);
+ //printf("Eagle cooling.h 1d H_plus_He_cooling electron abundance %.5e %.5e\n", temp_lambda, h_plus_he_electron_abundance);
if (element_lambda != NULL) element_lambda[0] = temp_lambda;
/* ------------------ */
@@ -1354,27 +1658,97 @@ eagle_metal_cooling_rate_1d_table(
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, temp_i, d_temp) *
- (h_plus_he_electron_abundance / solar_electron_abundance) *
- solar_ratio[i+2]; // hydrogen and helium aren't included here.
- cooling_rate += temp_lambda;
- *dlambda_du +=
- (((double)element_cooling_table[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[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;
- }
+ //if (element_lambda == NULL) {
+ // temp_lambda = interpol_1d_dbl(element_cooling_table, temp_i, d_temp) *
+ // (h_plus_he_electron_abundance / solar_electron_abundance); // hydrogen and helium aren't included here.
+ // *dlambda_du =
+ // (((double)element_cooling_table[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[temp_i]) *
+ // ((double)H_plus_He_electron_abundance_table[temp_i]) /
+ // ((double)element_electron_abundance_table[temp_i])) /
+ // ((double)du);
+ //} else {
+ 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;
+ //printf("Eagle cooling.h 1d element cooling electron abundance, solar_ratio %d %.5e %.5e %.5e\n", i, temp_lambda, solar_electron_abundance, solar_ratio[i+2]);
+ }
+ //}
//printf("Eagle cooling.h dlambda_du du %.5e %.5e\n", *dlambda_du, du);
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;//, lambda_net1 = 0.0, lambda_net2 = 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);
+
+ //int u_i;
+ //float du;
+ //get_index_1d(cooling->Therm, cooling->N_Temp, logu/log_10, &u_i, &du);
+ //float u1 = pow(10.0,cooling->Therm[u_i]);
+ //float u2 = pow(10.0,cooling->Therm[u_i+1]);
+ //lambda_net1 = eagle_metal_cooling_rate(
+ // log10(u1), z_index, dz, n_h_i, d_n_h,
+ // He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
+ //lambda_net2 = eagle_metal_cooling_rate(
+ // log10(u2), z_index, dz, n_h_i, d_n_h,
+ // He_i, d_He, p, cooling, cosmo, internal_const, element_lambda, abundance_ratio);
+ //*dLambdaNet_du = (lambda_net2 - lambda_net1)/(u2 - u1);
+
+ return lambda_net;
+}
+
/*
* @brief Wrapper function previously used to calculate cooling rate and
* dLambda_du
@@ -1423,6 +1797,66 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(
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
+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;
+ //const double log_10 = 2.30258509299404;
+ 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.
@@ -1448,7 +1882,7 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate_1d_table(
__attribute__((always_inline)) INLINE static double
eagle_print_metal_cooling_rate_1d_table(
double *H_plus_He_heat_table, double *H_plus_He_electron_abundance_table,
- double *element_cooling_table, double *element_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,
@@ -1456,15 +1890,16 @@ eagle_print_metal_cooling_rate_1d_table(
const struct phys_const *internal_const,
float *abundance_ratio) {
double *element_lambda, lambda_net = 0.0;
+ //const double log_10 = 2.30258509299404;
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[21];
+ 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_output_", element, ".dat");
+ 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");
@@ -1474,8 +1909,8 @@ eagle_print_metal_cooling_rate_1d_table(
for (int j = 0; j < cooling->N_Elements + 2; j++) element_lambda[j] = 0.0;
lambda_net = eagle_metal_cooling_rate_1d_table(
- log(u), &dLambdaNet_du, H_plus_He_heat_table,
- H_plus_He_electron_abundance_table, element_cooling_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++) {
@@ -1520,11 +1955,12 @@ __attribute__((always_inline)) INLINE static float bisection_iter(
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 *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;
- float shift = 0.3;
- double log10_e = 0.4342944819032518;
+ //float shift = 0.3;
+ //const double log10_e = 0.4342944819032518;
+ const float bisection_tolerance = 1.0e-8;
float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
@@ -1539,10 +1975,10 @@ __attribute__((always_inline)) INLINE static float bisection_iter(
/* set helium and hydrogen reheating term */
double LambdaTune = 0.0;
- if (cooling -> he_reion == 1) LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
+ if (cooling->he_reion == 1) LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
- x_a = x_init;
- x_b = cooling->Therm[0] / log10_e;
+ x_a = *ub;
+ x_b = *lb;
f_a = (LambdaTune / (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,
@@ -1554,17 +1990,17 @@ __attribute__((always_inline)) INLINE static float bisection_iter(
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);
- while (f_a * f_b >= 0 & i < 10 * eagle_max_iterations) {
- if ((x_a + shift) / log(10) < cooling->Therm[cooling->N_Temp - 1]) {
- x_a += shift;
- f_a = (LambdaTune / (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);
- }
- i++;
- }
+ //while (f_a * f_b >= 0 & i < 10 * eagle_max_iterations) {
+ // if ((x_a + shift) / log(10) < cooling->Therm[cooling->N_Temp - 1]) {
+ // x_a += shift;
+ // f_a = (LambdaTune / (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);
+ // }
+ // i++;
+ //}
#if SWIFT_DEBUG_CHECKS
assert(f_a * f_b < 0);
#endif
@@ -1585,190 +2021,13 @@ __attribute__((always_inline)) INLINE static float bisection_iter(
f_a = f_next;
}
i++;
- } while (fabs(x_a - x_b) > 1.0e-2 && i < 50*eagle_max_iterations);
+ } 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 Work in progress more stable integration scheme (brent-dekker)
- * that may be used for particles which do not
- * converge using newton_iter (may be faster than bisection_iter)
- *
- * @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 brent_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,
- 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) {
- /* this routine does the iteration scheme, call one and it iterates to
- * convergence */
-
- double x_a, x_b, x_c, x_d, x_s, x_temp, delta = 1.0e-2;
- double dLambdaNet_du, Lambda_a, Lambda_b, Lambda_c, Lambda_s;
- double f_a, f_b, f_c, f_s, f_temp;
- 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 */
- float LambdaTune = eagle_helium_reionization_extraheat(z_index, dz, cooling);
- static double zold = 100;
- if (zold > z_index) {
- float LambdaCumul = 0.0;
- LambdaCumul += LambdaTune;
- printf(" EagleDoCooling %d %g %g %g\n", z_index, dz, LambdaTune, LambdaCumul);
- zold = z_index;
- }
-
- int i = 0;
- x_a = x_init;
- // x_b = cooling->Therm[0]*log(10);
- x_b = 0.5 * x_init;
- x_c = x_a;
- x_d = 0.0;
- Lambda_a = 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);
- Lambda_b = 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);
- f_a = exp(x_a) - exp(x_init) - Lambda_a * ratefact * dt;
- f_b = exp(x_b) - exp(x_init) - Lambda_b * ratefact * dt;
- assert(f_a * f_b < 0);
- if (f_a < f_b) {
- x_temp = x_a;
- x_a = x_b;
- x_b = x_temp;
- f_temp = f_a;
- f_a = f_b;
- f_b = f_temp;
- }
- int mflag = 1;
-
- do /* iterate to convergence */
- {
- Lambda_c = eagle_cooling_rate_1d_table(
- x_c, &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_c = exp(x_c) - exp(x_init) - Lambda_c * ratefact * dt;
-
- if ((f_a != f_c) && (f_b != f_c)) {
- x_s = x_a * f_b * f_c / ((f_a - f_b) * (f_a - f_c)) +
- x_b * f_a * f_c / ((f_b - f_a) * (f_b - f_c)) +
- x_c * f_b * f_a / ((f_c - f_b) * (f_c - f_a));
- printf("Eagle cooling.h 1 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
- p->id, x_s, exp(x_s), x_a - x_b, i);
- printf(
- "Eagle cooling.h 1 id u_a, u_b, u_c, f_a, f_b, f_c, %llu %.5e %.5e "
- "%.5e %.5e %.5e %.5e %.5e %d\n",
- p->id, exp(x_a), exp(x_b), exp(x_c), f_a, f_b, f_c, x_a - x_b, i);
- } else {
- x_s = x_b - f_b * (x_b - x_a) / (f_b - f_a);
- printf("Eagle cooling.h 2 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
- p->id, x_s, exp(x_s), x_a - x_b, i);
- }
-
- if ((x_s < (3.0 * x_a + x_b) / 4.0 || x_s > x_b) ||
- (mflag == 1 && fabs(x_s - x_b) >= 0.5 * fabs(x_b - x_c)) ||
- (mflag != 1 && fabs(x_s - x_b) >= 0.5 * fabs(x_c - x_d)) ||
- (mflag == 1 && fabs(x_b - x_c) < delta) ||
- (mflag != 1 && fabs(x_c - x_d) < delta)) {
- x_s = 0.5 * (x_a + x_b);
- printf("Eagle cooling.h 3 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
- p->id, x_s, exp(x_s), x_a - x_b, i);
- mflag = 1;
- } else {
- mflag = 0;
- }
- x_d = x_c;
- x_c = x_b;
- printf("Eagle cooling.h 4 id x_s, u_s, error, %llu %.5e %.5e %.5e %d\n",
- p->id, x_s, exp(x_s), x_a - x_b, i);
- Lambda_s = eagle_cooling_rate_1d_table(
- x_s, &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_s = exp(x_s) - exp(x_init) - Lambda_s * ratefact * dt;
- if (f_s * f_a < 0) {
- x_b = x_s;
- } else {
- x_a = x_s;
- }
- Lambda_a = 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);
- Lambda_b = 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);
- f_a = exp(x_a) - exp(x_init) - Lambda_a * ratefact * dt;
- f_b = exp(x_b) - exp(x_init) - Lambda_b * ratefact * dt;
- if (f_a < f_b) {
- x_temp = x_a;
- x_a = x_b;
- x_b = x_temp;
- f_temp = f_a;
- f_a = f_b;
- f_b = f_temp;
- }
- i++;
-
- } while (fabs(x_b - x_a) > delta && i < eagle_max_iterations);
-
- if (i >= eagle_max_iterations) {
- n_eagle_cooling_rate_calls_3++;
- printf(
- "Eagle cooling.h particle %llu with density %.5e, initial energy, %.5e "
- "and redshift %.5e not converged \n",
- p->id, inn_h, exp(x_init), cosmo->z);
- }
-
- return x_s;
-}
/*
* @brief Newton Raphson integration scheme to calculate particle cooling over
@@ -1795,21 +2054,21 @@ __attribute__((always_inline)) INLINE static float brent_iter(
* @param phys_const Physical constants structure
* @param dt Hydro timestep
* @param printflag Flag to identify if scheme failed to converge
- * @param relax_factor Factor to increase stability
+ * @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, 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 x_init, double u_ini, int z_index,
float dz, int n_h_i, float d_n_h, int He_i, float d_He,
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 relax_factor) {
+ 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];
@@ -1836,51 +2095,44 @@ __attribute__((always_inline)) INLINE static float newton_iter(
x_old = x_init;
x = x_old;
int i = 0;
+ float relax = 1.0;//, relax_factor = 0.5;
+ float log_table_bound_high = 43.749, log_table_bound_low = 20.723;
do /* iterate to convergence */
{
- x_old = x;
+ if (relax == 1.0) x_old = x;
LambdaNet = (LambdaTune / (dt * ratefact)) +
- eagle_cooling_rate_1d_table(
- x_old, &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,
+ 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 iterate the variable.
+ // Newton iteration.
x = x_old -
- relax_factor * (1.0 - u_ini * exp(-x_old) -
+ relax * (1.0 - u_ini * exp(-x_old) -
LambdaNet * ratefact * dt * exp(-x_old)) /
(1.0 - dLambdaNet_du * ratefact * dt);
- // try limiting how far iterations can jump by
- if (*printflag == 1 && x < 0.3 * x_old) {
- x = 0.3 * x_old;
- } else if (*printflag == 1 && x > 1.3 * x_old) {
- x = 1.3 * x_old;
- }
// check whether iterations go out of bounds of table,
- // can be used to try to prevent runaway
- // Remove hardcoded bounds (log(1e19) = 43.749, log(1e9) = 20.723)!
- int table_out_of_bounds = 0;
- if (x > 43.749) {
- // x = log(u_ini);
- table_out_of_bounds = 1;
- }
- if (x < 20.723) {
- // x = log(u_ini);
- table_out_of_bounds = 1;
- }
-
- //residual = exp(x_old) - u_ini - LambdaNet * ratefact * dt;
- if (dt > 0 && *printflag == 1) {
- table_out_of_bounds = 0;
+ // 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;
}
+ // relax = relax_factor;
+ //} else {
+ // relax = 1.0;
+ //}
i++;
- if (table_out_of_bounds == 1) i = eagle_max_iterations;
- } while (fabs(x - x_old) > 1.0e-2 && i < eagle_max_iterations);
+ } 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
@@ -1918,6 +2170,80 @@ __attribute__((always_inline)) INLINE static void abundance_ratio_to_solar(
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);
+ //construct_1d_table_from_3d(p, cooling, cosmo, phys_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);
+
+ // 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.
*
@@ -1936,6 +2262,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
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) *
cooling->internal_energy_scale;
float dz;
@@ -1945,7 +2273,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
float XH, HeFrac;
double inn_h;
- double ratefact, u, LambdaNet, LambdaTune = 0, dLambdaNet_du;
+ double ratefact, u, LambdaNet, LambdaNet_eq, LambdaTune = 0, dLambdaNet_du;
static double zold = 100, LambdaCumul = 0;
dt *= units_cgs_conversion_factor(us, UNIT_CONV_TIME);
@@ -1953,8 +2281,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
u = u_old;
double u_ini = u_old, u_temp;
- float *abundance_ratio;
- abundance_ratio = malloc((chemistry_element_count + 2)*sizeof(float));
+ 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];
@@ -1981,86 +2308,82 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
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);
+
- double temp_table[176]; // WARNING sort out how it is declared/allocated
- construct_1d_table_from_4d(p, cooling, cosmo, phys_const,
- cooling->table.element_cooling.temperature,
- z_index, dz, He_i, d_He, n_h_i, d_n_h, temp_table);
-
- // To be used when amount of cooling is small, put back in when ready
- LambdaNet =
- eagle_metal_cooling_rate(log(u_ini), temp_table, z_index, dz, n_h_i, d_n_h,
- He_i, d_He, p, cooling, cosmo, phys_const, NULL, abundance_ratio);
- if (fabs(ratefact * LambdaNet * dt) < 0.05 * u_old) {
- /* cooling rate is small, take the explicit solution */
- u = u_old + ratefact * LambdaNet * dt;
+ 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 {
- // construct 1d table of cooling rates wrt temperature
- double H_plus_He_heat_table[176];
- double H_plus_He_electron_abundance_table[176];
- double element_cooling_table[176];
- double element_electron_abundance_table[176];
-
- // element cooling
- construct_1d_table_from_4d_H_He(
- p, cooling, cosmo, phys_const,
- cooling->table.element_cooling.H_plus_He_heating, z_index, dz, He_i,
- d_He, n_h_i, d_n_h, H_plus_He_heat_table, &u_temp, u_ini, dt);
- construct_1d_table_from_4d(
- p, cooling, cosmo, phys_const,
- cooling->table.element_cooling.H_plus_He_electron_abundance, z_index,
- dz, He_i, d_He, n_h_i, d_n_h, H_plus_He_electron_abundance_table);
- construct_1d_table_from_4d_elements(
- p, cooling, cosmo, phys_const,
- cooling->table.element_cooling.metal_heating, z_index, dz, n_h_i, d_n_h,
- element_cooling_table);
- construct_1d_table_from_3d(
- p, cooling, cosmo, phys_const,
- cooling->table.element_cooling.electron_abundance, z_index, dz, n_h_i,
- d_n_h, element_electron_abundance_table);
-
- n_eagle_cooling_rate_calls_2++;
-
- LambdaNet = eagle_cooling_rate_1d_table(
- log(u_ini), &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);
-
- double u_eq = 5.0e12;
-
- if (dt > 0) {
- float u_temp_guess = u_ini + LambdaNet * ratefact * dt;
- // printf("Eagle cooling.h id u_temp_guess, u_ini, LambdaNet, ratefact, dt
- // %llu %.5e %.5e %.5e %.5e %.5e \n", p->id,u_temp_guess, u_ini,
- // LambdaNet, ratefact, dt);
- // if ((LambdaNet < 0 && u_temp < u_temp_guess) ||
- // (LambdaNet >= 0 && u_temp > u_temp_guess))
- // u_temp = u_temp_guess;
- u_temp = u_temp_guess;
- if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
- u_temp = u_eq;
- if (isnan(u_temp)) u_temp = u_eq;
- double log_u_temp = log(u_temp);
-
- int printflag = 0;
- float relax_factor = 1.0;
- float x =
- newton_iter(log_u_temp, u_ini, 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, cosmo, cooling, phys_const,
- abundance_ratio, dt, &printflag, relax_factor);
- if (printflag == 1) {
- // newton method didn't work, so revert to bisection
- x =
- bisection_iter(log_u_temp,u_ini,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,cosmo,cooling,phys_const,abundance_ratio,dt);
- //n_eagle_cooling_rate_calls_4++;
- printflag = 2;
- }
- u = exp(x);
- }
- }
+ 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 = 1.0e12;//, u_max = pow(10.0,cooling->Therm[cooling->N_Temp-1]);
+ //do {
+ u_eq = u_ini;
+ LambdaNet_eq = LambdaTune/(dt*ratefact) + eagle_cooling_rate(
+ log(u_eq), &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
+ //} while (LambdaNet_eq > 0 && u_eq < u_max);
+
+ if (LambdaNet_eq < 0 && LambdaNet < 0) {
+ upper_bound = fmin(log(u_eq), log(u_ini));
+ } else if (LambdaNet_eq*LambdaNet < 0) {
+ upper_bound = fmax(log(u_eq), log(u_ini));
+ lower_bound = fmin(log(u_eq), log(u_ini));
+ } else {
+ upper_bound = fmax(log(u_eq), log(u_ini));
+ }
+
+ if (dt > 0) {
+ u_temp = u_ini + LambdaNet * ratefact * dt;
+ if ((LambdaNet < 0 && u_temp < u_eq) || (LambdaNet >= 0 && u_temp > u_eq))
+ u_temp = u_eq;
+ if (isnan(u_temp)) u_temp = u_eq;
+ double log_u_temp = log(u_temp);
+
+ int printflag = 0;
+ x = newton_iter(log_u_temp, u_ini, z_index, dz,
+ n_h_i, d_n_h, He_i, d_He, p, cosmo, cooling, phys_const,
+ abundance_ratio, dt, &printflag, &lower_bound, &upper_bound);
+ if (printflag == 1) {
+ //double Lambda_upper = LambdaTune/(dt*ratefact) + eagle_cooling_rate(
+ // upper_bound, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ // He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
+ //double Lambda_lower = LambdaTune/(dt*ratefact) + eagle_cooling_rate(
+ // lower_bound, &dLambdaNet_du, z_index, dz, n_h_i, d_n_h,
+ // He_i, d_He, p, cooling, cosmo, phys_const, abundance_ratio);
+ //printf("Eagle cooling.h id upper lower lambda_upper lambda_lower %llu %.5e %.5e %.5e %.5e\n", p->id, exp(upper_bound), exp(lower_bound), Lambda_upper, Lambda_lower);
+ // 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);
+ 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,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);
@@ -2103,13 +2426,43 @@ __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 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);
+ //const float internal_energy = hydro_get_physical_internal_energy(p,cosmo);
+ //if (fabs(xp->cooling_data.radiated_energy) > 1*internal_energy){
+ // float logu = log(internal_energy*cooling->internal_energy_scale);
+ // double dLambdaNet_du;
+ // int z_index, He_i, n_h_i;
+ // float dz, d_He, d_n_h;
+ // 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, phys_const) *
+ // cooling->number_density_scale;
+ // double ratefact = inn_h * (XH / eagle_proton_mass_cgs);
+
+ // 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);
+ //
+ // float abundance_ratio[chemistry_element_count + 2];
+ // abundance_ratio_to_solar(p, cooling, abundance_ratio);
+
+ // float Lambda_net = eagle_cooling_rate(logu, &dLambdaNet_du, z_index,
+ // dz, n_h_i, d_n_h, He_i, d_He,
+ // p, cooling, cosmo, phys_const, abundance_ratio);
+ // float tstep = fabs(internal_energy*cooling->internal_energy_scale/(Lambda_net*ratefact));
+ // //printf("Eagle cooling.h id tstep u lambda ratefact %llu %.5e %.5e %.5e %.5e\n", p->id, tstep, internal_energy*cooling->internal_energy_scale, Lambda_net, ratefact);
+ // return max(1.0e3*tstep,2.0e-9);
+ //} else {
+ // return FLT_MAX;
+ //}
+
return FLT_MAX;
}
diff --git a/src/cooling/EAGLE/eagle_cool_tables.h b/src/cooling/EAGLE/eagle_cool_tables.h
index aa1019b940e93461c929e146d91f9a916ba66df8..0d1805058baf2e8732654afe0581ffaaf4cef87b 100644
--- a/src/cooling/EAGLE/eagle_cool_tables.h
+++ b/src/cooling/EAGLE/eagle_cool_tables.h
@@ -538,8 +538,8 @@ inline 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(
- specs, j, k, cooling->N_Elements, cooling->N_Temp,
- cooling->N_nH); // Redshift invariant table!!!
+ specs, k, j, cooling->N_Elements, cooling->N_nH,
+ cooling->N_Temp); // Redshift invariant table!!!
cooling_table.metal_heating[cooling_index] =
-net_cooling_rate[table_index];
}
@@ -571,8 +571,8 @@ inline struct cooling_tables_redshift_invariant get_photodis_table(
table_index = row_major_index_3d(i, j, k, cooling->N_He,
cooling->N_Temp, cooling->N_nH);
cooling_index =
- row_major_index_3d(i, j, k, cooling->N_He, cooling->N_Temp,
- cooling->N_nH); // Redshift invariant table!!!
+ row_major_index_3d(k, i, j, cooling->N_nH, cooling->N_He,
+ cooling->N_Temp); // Redshift invariant table!!!
cooling_table.H_plus_He_heating[cooling_index] =
-he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
@@ -687,8 +687,6 @@ inline struct cooling_tables get_cooling_table(
cooling_index = row_major_index_4d(
z_index, specs, i, j, cooling->N_Redshifts, cooling->N_Elements,
cooling->N_nH, cooling->N_Temp);
- // cooling_index =
- // row_major_index_3d(specs,j,i,cooling->N_Elements,cooling->N_Temp,cooling->N_nH);
cooling_table.metal_heating[cooling_index] =
-net_cooling_rate[table_index];
}
@@ -722,8 +720,6 @@ inline struct cooling_tables get_cooling_table(
cooling_index =
row_major_index_4d(z_index, k, i, j, cooling->N_Redshifts,
cooling->N_nH, cooling->N_He, cooling->N_Temp);
- // cooling_index =
- // row_major_index_3d(i,j,k,cooling->N_He,cooling->N_Temp,cooling->N_nH);
cooling_table.H_plus_He_heating[cooling_index] =
-he_net_cooling_rate[table_index];
cooling_table.H_plus_He_electron_abundance[cooling_index] =
@@ -745,8 +741,6 @@ inline struct cooling_tables get_cooling_table(
table_index = row_major_index_2d(i, j, cooling->N_Temp, cooling->N_nH);
cooling_index = row_major_index_3d(z_index, j, i, cooling->N_Redshifts,
cooling->N_nH, cooling->N_Temp);
- // cooling_index =
- // row_major_index_2d(i,j,cooling->N_Temp,cooling->N_nH);
cooling_table.electron_abundance[cooling_index] =
electron_abundance[table_index];
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index fed7e5dd19fef380d2f3e06f3f3fe58e792d497f..36db326cf9624a77e3f2b4d06b5f67542f1edc9e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -312,7 +312,6 @@ __attribute__((always_inline)) INLINE static void
hydro_set_comoving_internal_energy_dt(struct part *restrict p, float du_dt) {
p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
- //if (p->id == 5643798559995) printf("Gadget2 particle id entropy entropy_dt du_dt %llu %.5e %.5e %.5e\n", p->id, p->entropy, p->entropy_dt, du_dt);
}
/**