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

testing cooling box with eagle tables

parent 4d3465ba
Branches
Tags
1 merge request!593Eagle cooling
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 136 additions and 1183 deletions
examples/CoolingBox/analytical_test.py
+ 58
− 61
View file @ a4428e45
import sys
import matplotlib
matplotlib.use("Agg")
from pylab import *
......@@ -31,18 +30,26 @@ import numpy as np
import h5py as h5
import sys
# function for interpolating 1d table of cooling rates
def interpol_lambda(u_list,cooling_rate,u):
if u < u_list[0]:
def interpol_lambda(temp_list,cooling_rate,temp):
#print(temp,temp_list[0],temp_list[len(temp_list)-1])
if temp < temp_list[0]:
return[cooling_rate[0]]
if u > u_list[len(u_list)-1]:
if temp > temp_list[len(temp_list)-1]:
return[cooling_rate[len(cooling_rate)-1]]
j = 0
while u_list[j+1] < u:
while temp_list[j+1] < temp:
j += 1
cooling = cooling_rate[j]*(u_list[j+1]-u)/(u_list[j+1]-u_list[j]) + cooling_rate[j+1]*(u-u_list[j])/(u_list[j+1]-u_list[j])
cooling = cooling_rate[j]*(temp_list[j+1]-temp)/(temp_list[j+1]-temp_list[j]) + cooling_rate[j+1]*(temp-temp_list[j])/(temp_list[j+1]-temp_list[j])
return cooling
def convert_u_to_temp_sol(u,rho):
k_b = 1.38E-16 #boltzmann
m_p = 1.67e-24 #proton mass
gamma = 5.0/3.0
pressure = u*rho*(gamma - 1.0)
temp = pressure/(k_b*rho/(0.59*m_p))
return temp
# File containing the total energy
stats_filename = "./energy.txt"
......@@ -67,86 +74,76 @@ unit_mass = units.attrs["Unit mass in cgs (U_M)"]
unit_length = units.attrs["Unit length in cgs (U_L)"]
unit_time = units.attrs["Unit time in cgs (U_t)"]
unit_vel = unit_length/unit_time
#hubble_param = 0.71
hubble_param = 1.0
unit_length = unit_length/hubble_param
unit_mass = unit_mass/hubble_param
yHe = 0.28
temp_0 = 1.0e7
rho = rho*unit_mass/(unit_length**3)
# Read snapshots
if len(sys.argv) >= 4:
nsnap = int(sys.argv[5])
else:
print("Need to specify number of snapshots, defaulting to 100.")
nsnap = 100
npart = 4096
nsnap = 40
npart = 32768
u_snapshots_cgs = zeros(nsnap)
u_part_snapshots_cgs = zeros((nsnap,npart))
t_snapshots_cgs = zeros(nsnap)
scale_factor = zeros(nsnap)
for i in range(nsnap):
snap = h5.File("coolingBox_%0.4d.hdf5"%i,'r')
u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*(unit_length**2)/(unit_time**2)
u_part_snapshots_cgs[i][:] = snap["/PartType0/InternalEnergy"][:]*unit_length**2/(unit_time**2)
t_snapshots_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
scale_factor[i] = snap["/Header"].attrs["Scale-factor"]
# calculate analytic solution. Since cooling rate is constant,
# temperature decrease in linear in time.
# read Lambda and temperatures from table
internal_energy = []
temperature = []
cooling_rate = []
file_in = open('cooling_output.dat', 'r')
for line in file_in:
data = line.split()
internal_energy.append(float(data[0]))
temperature.append(float(data[0]))
cooling_rate.append(-float(data[1]))
tfinal = t_snapshots_cgs[nsnap-1]
tfirst = t_snapshots_cgs[0]
nt = nsnap*10
dt = (tfinal-tfirst)/nt
tfinal = 3.3*t_snapshots_cgs[nsnap-1]
nt = 1e4
dt = tfinal/nt
t_sol = np.zeros(int(nt))
u_sol = np.zeros(int(nt))
temp_sol = np.zeros(int(nt))
lambda_sol = np.zeros(int(nt))
u = np.mean(u_part_snapshots_cgs[0,:])/scale_factor[0]**2
u_sol[0] = u
t_sol[0] = tfirst
# integrate explicit ODE
u = np.mean(u_part_snapshots_cgs[0,:])
temp_sol[0] = convert_u_to_temp_sol(u,rho)
#print(u,temp_sol[0])
for j in range(int(nt-1)):
t_sol[j+1] = t_sol[j] + dt
Lambda_net = interpol_lambda(internal_energy,cooling_rate,u_sol[j])
if j < 10:
print(u,Lambda_net)
if int(sys.argv[4]) == 1:
nH = 0.702*rho/(m_p)/scale_factor[0]**3
ratefact = nH*0.702/m_p
else:
nH = 0.752*rho/(m_p)/scale_factor[0]**3
ratefact = nH*0.752/m_p
u_next = u - Lambda_net*ratefact*dt
u_sol[j+1] = u_next
Lambda_net = interpol_lambda(temperature,cooling_rate,temp_sol[j])
#u_next = (u*m_p - Lambda_net*rho/(0.59*m_p)*dt)/m_p
nH = 0.75*rho/(m_p)
u_next = u - Lambda_net*nH**2/rho*dt
#print(u_next, u, Lambda_net,rho/(0.59*m_p),dt)
print(Lambda_net,rho,dt,nH)
temp_sol[j+1] = convert_u_to_temp_sol(u_next,rho)
lambda_sol[j] = Lambda_net
u = u_next
u_sol = u_sol*scale_factor[0]**2
# swift solution
mean_u = np.zeros(nsnap)
mean_temp = np.zeros(nsnap)
for j in range(nsnap):
mean_u[j] = np.mean(u_part_snapshots_cgs[j,:])
# plot and save results
log_nh = float(sys.argv[2])
redshift = float(sys.argv[1])
p = plt.figure(figsize = (6,6))
p1, = plt.loglog(t_sol,u_sol,linewidth = 0.5,color = 'k', marker = '*', markersize = 1.5,label = 'explicit ODE')
p2, = plt.loglog(t_snapshots_cgs,mean_u,linewidth = 0.5,color = 'r', marker = '*', markersize = 1.5,label = 'swift mean')
mean_temp[j] = convert_u_to_temp_sol(np.mean(u_part_snapshots_cgs[j,:]),rho)
p = figure()
p1, = plot(t_sol,temp_sol,linewidth = 0.5,color = 'k',label = 'analytical')
p2, = plot(t_snapshots_cgs,mean_temp,linewidth = 0.5,color = 'r',label = 'swift mean')
l = legend(handles = [p1,p2])
xlabel("${\\rm{Time~[s]}}$", labelpad=0, fontsize = 14)
ylabel("Internal energy ${\\rm{[erg \cdot g^{-1}]}}$", fontsize = 14)
if int(sys.argv[4]) == 1:
title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, solar metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_solar.png"
elif int(sys.argv[4]) == 0:
title('$n_h = 10^{' + "{}".format(log_nh) + '} \\rm{cm}^{-3}, z = ' + "{}".format(redshift) + '$, zero metallicity,\n relative error: ' + "{0:.4f}".format( (u_sol[int(nt)-1] - mean_u[nsnap-1])/(u_sol[int(nt)-1])), fontsize = 16)
name = "z_"+str(sys.argv[1])+"_nh_"+str(sys.argv[2])+"_pressure_"+str(sys.argv[3])+"_zero_metal.png"
savefig(name, dpi=200)
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])))
xlabel("${\\rm{Time~[s]}}$", labelpad=0)
ylabel("${\\rm{Temperature~[K]}}$")
savefig("energy.png", dpi=200)
p = figure()
p1, = loglog(temp_sol,lambda_sol,linewidth = 0.5,color = 'k',label = 'analytical')
xlabel("${\\rm{Temperature~[K]}}$")
ylabel("${\\rm{Cooling~rate}}$")
savefig("cooling.png", dpi=200)
examples/CoolingBox/makeIC.py
+ 3
− 3
View file @ a4428e45
......@@ -27,8 +27,8 @@ from numpy import *
# Parameters
periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec
rho = 3450447.97588 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
P = 4455445544.55 # Pressure in code units (at 10^6K)
rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
P = 4.5e8 # Pressure in code units (at 10^7K)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
......@@ -36,7 +36,7 @@ fileName = "coolingBox.hdf5"
#---------------------------------------------------
# Read id, position and h from glass
glass = h5py.File("glassCube_16.hdf5", "r")
glass = h5py.File("glassCube_32.hdf5", "r")
ids = glass["/PartType0/ParticleIDs"][:]
pos = glass["/PartType0/Coordinates"][:,:] * boxSize
h = glass["/PartType0/SmoothingLength"][:] * boxSize
......
examples/CoolingBox/run.sh
+ 3
− 3
View file @ a4428e45
......@@ -2,7 +2,7 @@
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e glassCube_16.hdf5 ]
if [ ! -e glassCube_32.hdf5 ]
then
echo "Fetching initial glass file for the cooling box example..."
./getGlass.sh
......@@ -21,7 +21,7 @@ then
fi
# Run SWIFT
../swift -s -c -C -t 16 -n 1000 coolingBox.yml
../swift -s -C -t 4 coolingBox.yml
# Check energy conservation and cooling rate
# python energy_plot.py
python energy_plot.py
src/cooling/EAGLE/cooling.h
+ 57
− 24
View file @ a4428e45
......@@ -41,7 +41,6 @@
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
// #include "physical_constants_cgs.h"
#include "units.h"
#include "eagle_cool_tables.h"
......@@ -111,6 +110,8 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable
*dx = 1;
} else {
*i = (int)floor(((float)x - table[0]) * dxm1);
//printf("Eagle cooling.h ntable, i, x, table[0], table[i], table[ntable-1], dxm1 %d %d %.5e %.5e %.5e %.5e %.5e\n", ntable, *i, x, table[0], table[*i], table[ntable-1], dxm1);
fflush(stdout);
*dx = ((float)x - table[*i]) * dxm1;
}
}
......@@ -129,7 +130,7 @@ __attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int
*dx = 1;
} else {
*i = (int)floor(((float)x - table[0]) * dxm1);
printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
//printf("i, x, dxm1: %d, %f, %f, %f\n", *i, (float) x, table[0], dxm1);
*dx = ((float)x - table[*i]) * dxm1;
}
}
......@@ -184,6 +185,8 @@ __attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int
result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
//printf("Eagle cooling.h interpol_2d dx dy table indices (0-3) nx ny %.5e %.5e %d %d %d %d %d %d\n", dx, dy, index[0], index[1], index[2], index[3], nx, ny);
return result;
}
......@@ -318,7 +321,8 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int
inline int set_cooling_SolarAbundances(const float *element_abundance,
double *cooling_element_abundance,
const struct cooling_function_data* restrict cooling) {
const struct cooling_function_data* restrict cooling,
const struct part* restrict p) {
int i, index;
int Silicon_SPH_Index = -1;
int Calcium_SPH_Index = -1;
......@@ -372,30 +376,56 @@ inline int set_cooling_SolarAbundances(const float *element_abundance,
if (strcmp(chemistry_get_element_name((enum chemistry_element) i), "Silicon") == 0) sili_index = i;
}
// Eagle way of identifying and assigning element abundance with strange workaround for calcium and sulphur
//for (i = 0; i < cooling->N_Elements; i++) {
// if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
// /* SPH does not track Calcium: use Si abundance */
// if (Silicon_SPH_Index == -1)
// cooling_element_abundance[i] = 0.0;
// else{
// cooling_element_abundance[i] =
// element_abundance[Silicon_SPH_Index] *
// cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
// }
// else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
// /* SPH does not track Sulphur: use Si abundance */
// if (Silicon_SPH_Index == -1)
// cooling_element_abundance[i] = 0.0;
// else{
// cooling_element_abundance[i] =
// element_abundance[Silicon_SPH_Index] *
// cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
// }
// else{
// cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
// cooling_ElementAbundance_SOLARM1[i];
// //printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
// }
//}
for (i = 0; i < cooling->N_Elements; i++) {
if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1)
/* SPH does not track Calcium: use Si abundance */
if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index != -1)
if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0;
else{
cooling_element_abundance[i] =
element_abundance[Silicon_SPH_Index] *
cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
element_abundance[Silicon_SPH_Index] * cooling->calcium_over_silicon_ratio *
cooling_ElementAbundance_SOLARM1[Calcium_CoolHeat_Index];
}
else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index == -1)
else if (i == Sulphur_CoolHeat_Index && Sulphur_SPH_Index != -1)
/* SPH does not track Sulphur: use Si abundance */
if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0;
else{
cooling_element_abundance[i] =
element_abundance[Silicon_SPH_Index] *
cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index];
element_abundance[Silicon_SPH_Index] * cooling->sulphur_over_silicon_ratio *
cooling_ElementAbundance_SOLARM1[Sulphur_CoolHeat_Index];
}
else{
cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
cooling_ElementAbundance_SOLARM1[i];
//printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i]);
}
//printf ("Eagle cooling.h element, name, abundance, solar abundance, solarm1, cooling abundance %d %s %.5e %.5e %.5e %.5e %.5e\n",cooling->ElementNamePointers[i],cooling->ElementNames[i], element_abundance[cooling->ElementNamePointers[i]],cooling->SolarAbundances[i], cooling_ElementAbundance_SOLARM1[i], cooling_element_abundance[i], element_abundance[i]);
}
return 0;
......@@ -475,10 +505,13 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons
get_index_1d(cooling->nH, cooling->N_nH, log10(inn_h), &n_h_i, &d_n_h);
get_index_1d(cooling->Therm, cooling->N_Temp, log10(u), &u_i, &d_u);
//logT = interpol_2d(cooling->table.collisional_cooling.temperature, u_i, He_i, d_u, d_He, cooling->N_Temp, cooling->N_He);
logT = interpol_4d(cooling->table.element_cooling.temperature, 0, He_i, n_h_i,
u_i, dz, d_He, d_n_h, d_u, cooling->N_Redshifts,cooling->N_He,cooling->N_nH,cooling->N_Temp);
T = pow(10.0, logT);
if (u_i == 0 && d_u == 0) T *= u / pow(10.0, cooling->Therm[0]);
......@@ -530,10 +563,10 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* ------------------ */
LambdaNet =
interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, 0, He_i, n_h_i,
interpol_4d(cooling->table.collisional_cooling.H_plus_He_heating, z_index, He_i, n_h_i,
temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
h_plus_he_electron_abundance =
interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, 0, He_i, n_h_i,
interpol_4d(cooling->table.collisional_cooling.H_plus_He_electron_abundance, z_index, He_i, n_h_i,
temp_i, dz, d_He, d_n_h, d_temp,1,cooling->N_He,cooling->N_nH,cooling->N_Temp);
/* ------------------ */
......@@ -558,15 +591,15 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* for each element, find the abundance and multiply it
by the interpolated heating-cooling */
set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling);
set_cooling_SolarAbundances(p->chemistry_data.metal_mass_fraction, cooling_solar_abundances, cooling, p);
solar_electron_abundance =
interpol_3d(cooling->table.element_cooling.electron_abundance, 0, n_h_i, temp_i, dz,
interpol_3d(cooling->table.element_cooling.electron_abundance, z_index, n_h_i, temp_i, dz,
d_n_h, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); /* ne/n_h */
double element_cooling_lambda;
for (i = 0; i < cooling->N_Elements; i++){
element_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, 0, i, n_h_i,
element_cooling_lambda = interpol_4d(cooling->table.element_cooling.metal_heating, z_index, i, n_h_i,
temp_i, dz, 0.0, d_n_h, d_temp,cooling->N_Redshifts,cooling->N_Elements,cooling->N_nH,cooling->N_Temp) *
(h_plus_he_electron_abundance / solar_electron_abundance) *
cooling_solar_abundances[i];
......@@ -610,6 +643,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
int i;
static double zold = 100, LambdaCumul = 0;
dt *= units_cgs_conversion_factor(us,UNIT_CONV_TIME);
u = u_old;
u_lower = u;
......@@ -620,6 +654,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* convert Hydrogen mass fraction in Hydrogen number density */
inn_h = chemistry_get_number_density(p,cosmo,chemistry_element_H,phys_const)*cooling->number_density_scale;
//inn_h = hydro_get_physical_density(p,cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY) * XH / eagle_proton_mass_cgs;
/* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */
ratefact = inn_h * (XH / eagle_proton_mass_cgs);
......@@ -707,16 +742,12 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
}
float cooling_du_dt = 0.0;
if (dt > 0) cooling_du_dt = (u - u_old)/dt/cooling->power_scale;
const float hydro_du_dt = hydro_get_internal_energy_dt(p);
/* Integrate cooling equation to enforce energy floor */
if (u_old + hydro_du_dt * dt < cooling->min_energy) {
cooling_du_dt = 0.f;
} else if (u_old + (hydro_du_dt + cooling_du_dt) * dt < cooling->min_energy) {
cooling_du_dt = (u_old + dt * hydro_du_dt - cooling->min_energy) / dt;
if (dt > 0){
cooling_du_dt = (u - u_old)/dt/cooling->power_scale;
}
const float hydro_du_dt = hydro_get_internal_energy_dt(p);
/* Update the internal energy time derivative */
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
......@@ -804,6 +835,8 @@ static INLINE void cooling_init_backend(
parser_get_param_string(parameter_file, "EagleCooling:filename",cooling->cooling_table_path);
cooling->reionisation_redshift = parser_get_param_float(parameter_file, "EagleCooling:reionisation_redshift");
cooling->calcium_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:CalciumOverSilicon");
cooling->sulphur_over_silicon_ratio = parser_get_param_float(parameter_file, "EAGLEChemistry:SulphurOverSilicon");
GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname,cooling);
......
src/cooling/EAGLE/cooling_struct.h
+ 2
− 0
View file @ a4428e45
......@@ -113,6 +113,8 @@ struct cooling_function_data {
double power_scale;
char cooling_table_path[500];
float reionisation_redshift;
float calcium_over_silicon_ratio;
float sulphur_over_silicon_ratio;
};
......
src/cooling/EAGLE/eagle_cool_tables.h
+ 1
− 9
View file @ a4428e45
......@@ -261,24 +261,15 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
net_cooling_rate);
status = H5Dclose(dataset);
//for(i = 0; i < cooling->N_Temp*cooling->N_nH; i++) printf("eagle_cool_tables.h i, net_cooling_rate %d, %.5e\n",i,net_cooling_rate[i]);
for (j = 0; j < cooling->N_Temp; j++){
for (k = 0; k < cooling->N_nH; k++){
table_index = row_major_index_2d(j,k,cooling->N_Temp,cooling->N_nH);
cooling_index = row_major_index_4d(0,specs,k,j,1,cooling->N_Elements,cooling->N_nH,cooling->N_Temp); //Redshift invariant table!!!
cooling_table.metal_heating[cooling_index] = -net_cooling_rate[table_index];
//printf("eagle_cool_tables.h j,k,j*cooling->N_nH+k,table_index,cooling_index,net cooling rate, cooling_table value %d, %d, %d, %d, %d %.5e, %.5e\n",j,k,j*cooling->N_nH+k,table_index,cooling_index,-net_cooling_rate[table_index],cooling_table.metal_heating[cooling_index]);
}
}
}
//for (i = 0; i < cooling->N_nH; i++){
// table_index = row_major_index_2d(0,i,cooling->N_Temp,cooling->N_nH);
// cooling_index = row_major_index_4d(0,0,i,0,1,cooling->N_Elements,cooling->N_nH,cooling->N_Temp);
// printf("eagle_cool_tables.h i, cooling_index, table_index, cooling_table.metal_heating, net heating = %d %d %d %.5e %.5e\n",i,cooling_index,table_index,cooling_table.metal_heating[cooling_index],net_cooling_rate[table_index]);
//}
//for (i = 0; i < cooling->N_Elements*cooling->N_Temp*cooling->N_nH; i++) printf("eagle cooling.h i, cooling_table %d, %.5e\n",i,cooling_table.metal_heating[i]);
//error("stop in eagle_cool_tables.h");
/* Helium */
......@@ -382,6 +373,7 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool
cooling_table.H_plus_He_electron_abundance = (float *)malloc(cooling->N_He*cooling->N_Temp*cooling->N_nH*sizeof(float));
sprintf(fname, "%sz_photodis.hdf5", cooling_table_path);
//sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
......
src/cooling/EAGLE/eagle_init_cool.h deleted 100644 → 0
+ 0
− 899
View file @ 4d3465ba
This diff is collapsed.
src/cooling/EAGLE/eagle_interpol.h deleted 100644 → 0
+ 0
− 179
View file @ 4d3465ba
#ifndef SWIFT_EAGLE_INTERPOL_H
#define SWIFT_EAGLE_INTERPOL_H
#include "cooling_read_table.h"
/*
* ----------------------------------------------------------------------
* This routine returns the position i of a value x in a 1D table and the
* displacement dx needed for the interpolation. The table is assumed to
* be evenly spaced.
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable, double x, int *i, float *dx) {
float dxm1;
const float EPS = 1.e-4;
dxm1 = (float)(ntable - 1) / (table[ntable - 1] - table[0]);
if ((float)x <= table[0] + EPS) {
*i = 0;
*dx = 0;
} else if ((float)x >= table[ntable - 1] - EPS) {
*i = ntable - 2;
*dx = 1;
} else {
*i = (int)floor(((float)x - table[0]) * dxm1);
*dx = ((float)x - table[*i]) * dxm1;
}
}
/*
* ----------------------------------------------------------------------
* This routine performs a linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE float interpol_1d(float *table, int i, float dx) {
float result;
result = (1 - dx) * table[i] + dx * table[i + 1];
return result;
}
/*
* ----------------------------------------------------------------------
* This routine performs a linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE double interpol_1d_dbl(double *table, int i, float dx) {
double result;
result = (1 - dx) * table[i] + dx * table[i + 1];
return result;
}
/*
* ----------------------------------------------------------------------
* This routine performs a bi-linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int j, float dx, float dy, int size) {
float result;
int index[4];
index[0] = row_major_index_2d(i,j,size);
index[1] = row_major_index_2d(i,j+1,size);
index[2] = row_major_index_2d(i+1,j,size);
index[3] = row_major_index_2d(i+1,j+1,size);
result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
return result;
}
/*
* ----------------------------------------------------------------------
* This routine performs a bi-linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE double interpol_2d_dbl(double *table, int i, int j, double dx, double dy, int size) {
double result;
int index[4];
index[0] = row_major_index_2d(i,j,size);
index[1] = row_major_index_2d(i,j+1,size);
index[2] = row_major_index_2d(i+1,j,size);
index[3] = row_major_index_2d(i+1,j+1,size);
result = (1 - dx) * (1 - dy) * table[index[0]] + (1 - dx) * dy * table[index[1]] +
dx * (1 - dy) * table[index[2]] + dx * dy * table[index[3]];
return result;
}
/*
* ----------------------------------------------------------------------
* This routine performs a tri-linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE float interpol_3d(float *table, int i, int j, int k, float dx, float dy,
float dz, int size) {
float result;
int index[8];
index[0] = row_major_index_3d(i,j,k,size);
index[1] = row_major_index_3d(i,j,k+1,size);
index[2] = row_major_index_3d(i,j+1,k,size);
index[3] = row_major_index_3d(i,j+1,k+1,size);
index[4] = row_major_index_3d(i+1,j,k,size);
index[5] = row_major_index_3d(i+1,j,k+1,size);
index[6] = row_major_index_3d(i+1,j+1,k,size);
index[7] = row_major_index_3d(i+1,j+1,k+1,size);
result = (1 - dx) * (1 - dy) * (1 - dz) * table[index[0]] +
(1 - dx) * (1 - dy) * dz * table[index[1]] +
(1 - dx) * dy * (1 - dz) * table[index[2]] +
(1 - dx) * dy * dz * table[index[3]] +
dx * (1 - dy) * (1 - dz) * table[index[4]] +
dx * (1 - dy) * dz * table[index[5]] +
dx * dy * (1 - dz) * table[index[6]] +
dx * dy * dz * table[index[7]];
return result;
}
/*
* ----------------------------------------------------------------------
* This routine performs a quadri-linear interpolation
* ----------------------------------------------------------------------
*/
__attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int j, int k, int l, float dx,
float dy, float dz, float dw, int size) {
float result;
int index[16];
index[0] = row_major_index_4d(i,j,k,l,size);
index[1] = row_major_index_4d(i,j,k,l+1,size);
index[2] = row_major_index_4d(i,j,k+1,l,size);
index[3] = row_major_index_4d(i,j,k+1,l+1,size);
index[4] = row_major_index_4d(i,j+1,k,l,size);
index[5] = row_major_index_4d(i,j+1,k,l+1,size);
index[6] = row_major_index_4d(i,j+1,k+1,l,size);
index[7] = row_major_index_4d(i,j+1,k+1,l+1,size);
index[8] = row_major_index_4d(i+1,j,k,l,size);
index[9] = row_major_index_4d(i+1,j,k,l+1,size);
index[10] = row_major_index_4d(i+1,j,k+1,l,size);
index[11] = row_major_index_4d(i+1,j,k+1,l+1,size);
index[12] = row_major_index_4d(i+1,j+1,k,l,size);
index[13] = row_major_index_4d(i+1,j+1,k,l+1,size);
index[14] = row_major_index_4d(i+1,j+1,k+1,l,size);
index[15] = row_major_index_4d(i+1,j+1,k+1,l+1,size);
result = (1 - dx) * (1 - dy) * (1 - dz) * (1 - dw) * table[index[0]] +
(1 - dx) * (1 - dy) * (1 - dz) * dw * table[index[1]] +
(1 - dx) * (1 - dy) * dz * (1 - dw) * table[index[2]] +
(1 - dx) * (1 - dy) * dz * dw * table[index[3]] +
(1 - dx) * dy * (1 - dz) * (1 - dw) * table[index[4]] +
(1 - dx) * dy * (1 - dz) * dw * table[index[5]] +
(1 - dx) * dy * dz * (1 - dw) * table[index[6]] +
(1 - dx) * dy * dz * dw * table[index[7]] +
dx * (1 - dy) * (1 - dz) * (1 - dw) * table[index[8]] +
dx * (1 - dy) * (1 - dz) * dw * table[index[9]] +
dx * (1 - dy) * dz * (1 - dw) * table[index[10]] +
dx * (1 - dy) * dz * dw * table[index[11]] +
dx * dy * (1 - dz) * (1 - dw) * table[index[12]] +
dx * dy * (1 - dz) * dw * table[index[13]] +
dx * dy * dz * (1 - dw) * table[index[14]] +
dx * dy * dz * dw * table[index[15]];
return result;
}
#endif
tests/testCooling.c
+ 12
− 5
View file @ a4428e45
......@@ -68,15 +68,24 @@ int main() {
exit(1);
}
gamma = 5.0/3.0;
chemistry_init(params, &us, &internal_const, &chemistry_data);
chemistry_first_init_part(&p,&xp,&chemistry_data);
chemistry_print(&chemistry_data);
u = 1.0*pow(10.0,11)/(units_cgs_conversion_factor(&us,UNIT_CONV_ENERGY)/units_cgs_conversion_factor(&us,UNIT_CONV_MASS));
pressure = u*p.rho*(gamma -1.0);
hydrogen_number_density = hydrogen_number_density_cgs*pow(units_cgs_conversion_factor(&us,UNIT_CONV_LENGTH),3);
p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
p.entropy = pressure/(pow(p.rho,gamma));
cosmology_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo);
cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&cooling);
chemistry_init(params, &us, &internal_const, &chemistry_data);
chemistry_first_init_part(&p,&xp,&chemistry_data);
chemistry_print(&chemistry_data);
for(int t_i = 0; t_i < n_t_i; t_i++){
......@@ -86,8 +95,6 @@ int main() {
p.rho = hydrogen_number_density*internal_const.const_proton_mass*(1.0+p.chemistry_data.metal_mass_fraction[EAGLE_Helium]/p.chemistry_data.metal_mass_fraction[EAGLE_Hydrogen]);
p.entropy = pressure/(pow(p.rho,gamma));
gamma = 5.0/3.0;
cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const);
temperature_cgs = eagle_convert_u_to_temp(&p,&cooling,&cosmo,&internal_const);
fprintf(output_file,"%.5e %.5e\n",temperature_cgs,cooling_du_dt);
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment