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

testing cooling box with eagle tables

parent 261af9c3
Branches
Tags
1 merge request!593Eagle cooling
Hide whitespace changes
Inline Side-by-side
Showing
with 225 additions and 1120 deletions
examples/CoolingBox/analytical_test.py 0 → 100644
+ 149
− 0
View file @ 89c22621
import matplotlib
matplotlib.use("Agg")
from pylab import *
import h5py
# Plot parameters
params = {'axes.labelsize': 10,
'axes.titlesize': 10,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 10,
'ytick.labelsize': 10,
'text.usetex': True,
'figure.figsize' : (3.15,3.15),
'figure.subplot.left' : 0.2,
'figure.subplot.right' : 0.99,
'figure.subplot.bottom' : 0.2,
'figure.subplot.top' : 0.9,
'figure.subplot.wspace' : 0.15,
'figure.subplot.hspace' : 0.12,
'lines.markersize' : 6,
'lines.linewidth' : 3.,
'text.latex.unicode': True
}
rcParams.update(params)
rc('font',**{'family':'sans-serif','sans-serif':['Times']})
import numpy as np
import h5py as h5
import sys
def interpol_lambda(temp_list,cooling_rate,temp):
#print(temp,temp_list[0],temp_list[len(temp_list)-1])
if temp < temp_list[0]:
return[cooling_rate[0]]
if temp > temp_list[len(temp_list)-1]:
return[cooling_rate[len(cooling_rate)-1]]
j = 0
while temp_list[j+1] < temp:
j += 1
cooling = cooling_rate[j]*(temp_list[j+1]-temp)/(temp_list[j+1]-temp_list[j]) + cooling_rate[j+1]*(temp-temp_list[j])/(temp_list[j+1]-temp_list[j])
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"
# First snapshot
snap_filename = "coolingBox_0000.hdf5"
# Some constants in cgs units
k_b = 1.38E-16 #boltzmann
m_p = 1.67e-24 #proton mass
# Initial conditions set in makeIC.py
T_init = 1.0e7
# Read the initial state of the gas
f = h5.File(snap_filename,'r')
rho = np.mean(f["/PartType0/Density"])
pressure = np.mean(f["/PartType0/Pressure"])
# Read the units parameters from the snapshot
units = f["InternalCodeUnits"]
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
nsnap = 40
npart = 32768
u_snapshots_cgs = zeros(nsnap)
u_part_snapshots_cgs = zeros((nsnap,npart))
t_snapshots_cgs = 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)
t_snapshots_cgs[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
temperature = []
cooling_rate = []
file_in = open('cooling_output.dat', 'r')
for line in file_in:
data = line.split()
temperature.append(float(data[0]))
cooling_rate.append(-float(data[1]))
tfinal = 3.3*t_snapshots_cgs[nsnap-1]
nt = 1e4
dt = tfinal/nt
t_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,:])
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(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
mean_temp = 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')
l = legend(handles = [p1,p2])
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/coolingBox.yml
+ 2
− 2
View file @ 89c22621
...@@ -17,7 +17,7 @@ Cosmology: ...@@ -17,7 +17,7 @@ Cosmology:
# Parameters governing the time integration # Parameters governing the time integration
TimeIntegration: TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units). time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 0.25 # The end time of the simulation (in internal units). time_end: 0.02 # The end time of the simulation (in internal units).
dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units). dt_min: 1e-5 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units). dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
...@@ -25,7 +25,7 @@ TimeIntegration: ...@@ -25,7 +25,7 @@ TimeIntegration:
Snapshots: Snapshots:
basename: coolingBox # Common part of the name of output files basename: coolingBox # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units) time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-2 # Time difference between consecutive outputs (in internal units) delta_time: 5e-4 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics # Parameters governing the conserved quantities statistics
Statistics: Statistics:
......
examples/CoolingBox/makeIC.py
+ 1
− 1
View file @ 89c22621
...@@ -28,7 +28,7 @@ from numpy import * ...@@ -28,7 +28,7 @@ from numpy import *
periodic= 1 # 1 For periodic box periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec boxSize = 1 # 1 kiloparsec
rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3) rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
P = 4.5e6 # Pressure in code units (at 10^5K) P = 4.5e8 # Pressure in code units (at 10^7K)
gamma = 5./3. # Gas adiabatic index gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5" fileName = "coolingBox.hdf5"
......
examples/CoolingBox/run.sh
+ 1
− 1
View file @ 89c22621
...@@ -21,7 +21,7 @@ then ...@@ -21,7 +21,7 @@ then
fi fi
# Run SWIFT # Run SWIFT
../swift -s -C -t 1 coolingBox.yml ../swift -s -C -t 4 coolingBox.yml
# Check energy conservation and cooling rate # Check energy conservation and cooling rate
python energy_plot.py python energy_plot.py
src/cooling/EAGLE/cooling.h
+ 57
− 24
View file @ 89c22621
...@@ -41,7 +41,6 @@ ...@@ -41,7 +41,6 @@
#include "parser.h" #include "parser.h"
#include "part.h" #include "part.h"
#include "physical_constants.h" #include "physical_constants.h"
// #include "physical_constants_cgs.h"
#include "units.h" #include "units.h"
#include "eagle_cool_tables.h" #include "eagle_cool_tables.h"
...@@ -111,6 +110,8 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable ...@@ -111,6 +110,8 @@ __attribute__((always_inline)) INLINE void get_index_1d(float *table, int ntable
*dx = 1; *dx = 1;
} else { } else {
*i = (int)floor(((float)x - table[0]) * dxm1); *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; *dx = ((float)x - table[*i]) * dxm1;
} }
} }
...@@ -129,7 +130,7 @@ __attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int ...@@ -129,7 +130,7 @@ __attribute__((always_inline)) INLINE void get_index_1d_therm(float *table, int
*dx = 1; *dx = 1;
} else { } else {
*i = (int)floor(((float)x - table[0]) * dxm1); *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; *dx = ((float)x - table[*i]) * dxm1;
} }
} }
...@@ -184,6 +185,8 @@ __attribute__((always_inline)) INLINE float interpol_2d(float *table, int i, int ...@@ -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]] + 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]]; 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; return result;
} }
...@@ -318,7 +321,8 @@ __attribute__((always_inline)) INLINE float interpol_4d(float *table, int i, int ...@@ -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, inline int set_cooling_SolarAbundances(const float *element_abundance,
double *cooling_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 i, index;
int Silicon_SPH_Index = -1; int Silicon_SPH_Index = -1;
int Calcium_SPH_Index = -1; int Calcium_SPH_Index = -1;
...@@ -372,30 +376,56 @@ inline int set_cooling_SolarAbundances(const float *element_abundance, ...@@ -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; 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++) { for (i = 0; i < cooling->N_Elements; i++) {
if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index == -1) if (i == Calcium_CoolHeat_Index && Calcium_SPH_Index != -1)
/* SPH does not track Calcium: use Si abundance */
if (Silicon_SPH_Index == -1) if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0; cooling_element_abundance[i] = 0.0;
else{ else{
cooling_element_abundance[i] = cooling_element_abundance[i] =
element_abundance[Silicon_SPH_Index] * element_abundance[Silicon_SPH_Index] * cooling->calcium_over_silicon_ratio *
cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index]; 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 */ /* SPH does not track Sulphur: use Si abundance */
if (Silicon_SPH_Index == -1) if (Silicon_SPH_Index == -1)
cooling_element_abundance[i] = 0.0; cooling_element_abundance[i] = 0.0;
else{ else{
cooling_element_abundance[i] = cooling_element_abundance[i] =
element_abundance[Silicon_SPH_Index] * element_abundance[Silicon_SPH_Index] * cooling->sulphur_over_silicon_ratio *
cooling_ElementAbundance_SOLARM1[Silicon_CoolHeat_Index]; cooling_ElementAbundance_SOLARM1[Sulphur_CoolHeat_Index];
} }
else{ else{
cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] * cooling_element_abundance[i] = element_abundance[cooling->ElementNamePointers[i]] *
cooling_ElementAbundance_SOLARM1[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; return 0;
...@@ -475,10 +505,13 @@ __attribute__((always_inline)) INLINE static double eagle_convert_u_to_temp(cons ...@@ -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->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); 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, 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); 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); 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) T *= u / pow(10.0, cooling->Therm[0]);
...@@ -530,10 +563,10 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str ...@@ -530,10 +563,10 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* ------------------ */ /* ------------------ */
LambdaNet = 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); 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 = 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); 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 ...@@ -558,15 +591,15 @@ __attribute__((always_inline)) INLINE static double eagle_cooling_rate(const str
/* for each element, find the abundance and multiply it /* for each element, find the abundance and multiply it
by the interpolated heating-cooling */ 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 = 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 */ d_n_h, d_temp,cooling->N_Redshifts,cooling->N_nH,cooling->N_Temp); /* ne/n_h */
double element_cooling_lambda; double element_cooling_lambda;
for (i = 0; i < cooling->N_Elements; i++){ 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) * 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) *
cooling_solar_abundances[i]; cooling_solar_abundances[i];
...@@ -610,6 +643,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -610,6 +643,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
int i; int i;
static double zold = 100, LambdaCumul = 0; static double zold = 100, LambdaCumul = 0;
dt *= units_cgs_conversion_factor(us,UNIT_CONV_TIME);
u = u_old; u = u_old;
u_lower = u; u_lower = u;
...@@ -620,6 +654,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -620,6 +654,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* convert Hydrogen mass fraction in Hydrogen number density */ /* 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 = 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 /* ratefact = inn_h * inn_h / rho; Might lead to round-off error: replaced by
* equivalent expression below */ * equivalent expression below */
ratefact = inn_h * (XH / eagle_proton_mass_cgs); ratefact = inn_h * (XH / eagle_proton_mass_cgs);
...@@ -707,16 +742,12 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -707,16 +742,12 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
} }
float cooling_du_dt = 0.0; float cooling_du_dt = 0.0;
if (dt > 0) cooling_du_dt = (u - u_old)/dt/cooling->power_scale; if (dt > 0){
const float hydro_du_dt = hydro_get_internal_energy_dt(p); cooling_du_dt = (u - u_old)/dt/cooling->power_scale;
/* 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;
} }
const float hydro_du_dt = hydro_get_internal_energy_dt(p);
/* Update the internal energy time derivative */ /* Update the internal energy time derivative */
hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt); hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
...@@ -804,6 +835,8 @@ static INLINE void cooling_init_backend( ...@@ -804,6 +835,8 @@ static INLINE void cooling_init_backend(
parser_get_param_string(parameter_file, "EagleCooling:filename",cooling->cooling_table_path); 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->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); GetCoolingRedshifts(cooling);
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path); sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
ReadCoolingHeader(fname,cooling); ReadCoolingHeader(fname,cooling);
......
src/cooling/EAGLE/cooling_struct.h
+ 2
− 0
View file @ 89c22621
...@@ -113,6 +113,8 @@ struct cooling_function_data { ...@@ -113,6 +113,8 @@ struct cooling_function_data {
double power_scale; double power_scale;
char cooling_table_path[500]; char cooling_table_path[500];
float reionisation_redshift; float reionisation_redshift;
float calcium_over_silicon_ratio;
float sulphur_over_silicon_ratio;
}; };
......
src/cooling/EAGLE/eagle_cool_tables.h
+ 1
− 9
View file @ 89c22621
...@@ -261,24 +261,15 @@ inline struct cooling_tables_redshift_invariant get_no_compt_table(char *cooling ...@@ -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, status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
net_cooling_rate); net_cooling_rate);
status = H5Dclose(dataset); 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 (j = 0; j < cooling->N_Temp; j++){
for (k = 0; k < cooling->N_nH; k++){ for (k = 0; k < cooling->N_nH; k++){
table_index = row_major_index_2d(j,k,cooling->N_Temp,cooling->N_nH); 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_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]; 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 */ /* Helium */
...@@ -382,6 +373,7 @@ inline struct cooling_tables_redshift_invariant get_collisional_table(char *cool ...@@ -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)); 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_photodis.hdf5", cooling_table_path);
//sprintf(fname, "%sz_collis.hdf5", cooling_table_path);
file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT); file_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
......
src/cooling/EAGLE/eagle_init_cool.h deleted 100644 → 0
+ 0
− 899
View file @ 261af9c3
#ifndef SWIFT_EAGLE_INIT_COOL_H
#define SWIFT_EAGLE_INIT_COOL_H
/*
* ----------------------------------------------------------------------
* Prototype for built-in testing routines
* ----------------------------------------------------------------------
*/
void DriverEagleTestCool();
void EagleTestCool(float z);
void DriverEagleTestCoolThermalEvolution();
/*
* ----------------------------------------------------------------------
* This routine gets the redshifts from a text file
* ----------------------------------------------------------------------
*/
void GetCoolingRedshifts() {
FILE *infile;
int i = 0;
char buffer[500], redfilename[500];
sprintf(redfilename, "%s/redshifts.dat", All.CoolTablePath);
infile = fopen(redfilename, "r");
if (infile == NULL) puts("GetCoolingRedshifts can't open a file");
if (fscanf(infile, "%s", buffer) != EOF) {
cooling_N_Redshifts = atoi(buffer);
cooling_Redshifts =
(float *)mymalloc("cNz", cooling_N_Redshifts * sizeof(float));
while (fscanf(infile, "%s", buffer) != EOF) {
cooling_Redshifts[i] = atof(buffer);
i += 1;
}
}
fclose(infile);
#ifdef EAGLE_VERBOSE
for (i = 0; i < cooling_N_Redshifts; i++)
printf("Cooling table %2d at redshift %1.3f\n", i + 1,
cooling_Redshifts[i]);
printf("\n");
#endif
}
/*
* ----------------------------------------------------------------------
* This routine reads in the header of the cooling table files
* ----------------------------------------------------------------------
*/
void ReadCoolingHeader(char *fname) {
int i;
hid_t tempfile_id, dataset, datatype;
herr_t status;
/* fill the constants */
tempfile_id = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
if (tempfile_id < 0) {
printf("[ReadCoolingHeader()]: unable to open file %s\n", fname);
endrun(101);
}
dataset = H5Dopen(tempfile_id, "/Header/Number_of_density_bins");
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling_N_nH);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Header/Number_of_temperature_bins");
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling_N_Temp);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Header/Number_of_helium_fractions");
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling_N_He);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Header/Number_of_metals");
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling_N_Elements);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Header/Abundances/Number_of_abundances");
status = H5Dread(dataset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
&cooling_N_SolarAbundances);
status = H5Dclose(dataset);
/* allocate arrays for cooling table header */
allocate_header_arrays();
/* fill the arrays */
dataset = H5Dopen(tempfile_id, "/Solar/Temperature_bins");
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_Temp);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Solar/Hydrogen_density_bins");
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_nH);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Metal_free/Temperature/Energy_density_bins");
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_Therm);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Metal_free/Helium_mass_fraction_bins");
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_HeFrac);
status = H5Dclose(dataset);
dataset = H5Dopen(tempfile_id, "/Header/Abundances/Solar_mass_fractions");
status = H5Dread(dataset, H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL, H5P_DEFAULT,
cooling_SolarAbundances);
status = H5Dclose(dataset);
char element_names[cooling_N_Elements][EL_NAME_LENGTH];
hsize_t string_length = EL_NAME_LENGTH;
/* names of chemical elements stored in table */
datatype = H5Tcopy(H5T_C_S1);
status = H5Tset_size(datatype, string_length);
dataset = H5Dopen(tempfile_id, "/Header/Metal_names");
status =
H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, element_names);
status = H5Dclose(dataset);
for (i = 0; i < cooling_N_Elements; i++)
cooling_ElementNames[i] = mystrdup(element_names[i]);
char solar_abund_names[cooling_N_SolarAbundances][EL_NAME_LENGTH];
/* assumed solar abundances used in constructing the tables, and corresponding
* names */
dataset = H5Dopen(tempfile_id, "/Header/Abundances/Abund_names");
status = H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT,
solar_abund_names);
status = H5Dclose(dataset);
H5Tclose(datatype);
for (i = 0; i < cooling_N_SolarAbundances; i++)
cooling_SolarAbundanceNames[i] = mystrdup(solar_abund_names[i]);
status = H5Fclose(tempfile_id);
/* Convert to temperature, density and internal energy arrays to log10 */
for (i = 0; i < cooling_N_Temp; i++) {
cooling_Temp[i] = log10(cooling_Temp[i]);
cooling_Therm[i] = log10(cooling_Therm[i]);
}
for (i = 0; i < cooling_N_nH; i++) cooling_nH[i] = log10(cooling_nH[i]);
printf("Done with cooling table header.\n");
fflush(stdout);
#ifdef EAGLE_VERBOSE
printf(
"\n\nTemperature grid (LOG10) for energy to temperature conversion\n\n");
for (i = 0; i < cooling_N_Temp; i++)
printf("Temperature[%3d] = %f\n", i, cooling_Temp[i]);
printf("\n\nEnergy grid (LOG10) for energy to temperature conversion\n\n");
for (i = 0; i < cooling_N_Temp; i++)
printf("Energy[%3d] = %f\n", i, cooling_Therm[i]);
printf("\n\nDensity grid (LOG10) for energy to temperature conversion\n\n");
for (i = 0; i < cooling_N_nH; i++)
printf("Density[%2d] = %f\n", i, cooling_nH[i]);
#endif
}
/*
* ----------------------------------------------------------------------
* This routine broadcasts the header infos
* ----------------------------------------------------------------------
*/
void BroadCastHeader() {
int position, bufsize, size, i;
char *buffer;
/* get size of the buffer */
MPI_Pack_size(6, MPI_INT, MPI_COMM_WORLD, &size);
bufsize = size;
/* allocate memory for the buffer */
buffer = (char *)mymalloc("buffer", bufsize);
if (ThisTask == 0) {
position = 0;
/* pack the array dimensions */
MPI_Pack(&cooling_N_nH, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(&cooling_N_Temp, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(&cooling_N_He, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(&cooling_N_Elements, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(&cooling_N_SolarAbundances, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(&cooling_N_Redshifts, 1, MPI_INT, buffer, bufsize, &position,
MPI_COMM_WORLD);
}
MPI_Bcast(buffer, bufsize, MPI_PACKED, 0, MPI_COMM_WORLD);
if (ThisTask != 0) {
position = 0;
/* unpack the array dimensions */
MPI_Unpack(buffer, bufsize, &position, &cooling_N_nH, 1, MPI_INT,
MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, &cooling_N_Temp, 1, MPI_INT,
MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, &cooling_N_He, 1, MPI_INT,
MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, &cooling_N_Elements, 1, MPI_INT,
MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, &cooling_N_SolarAbundances, 1,
MPI_INT, MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, &cooling_N_Redshifts, 1, MPI_INT,
MPI_COMM_WORLD);
}
/* free allocated memory */
myfree(buffer);
if (ThisTask != 0) {
/* allocate header arrays */
allocate_header_arrays();
/* allocate redshift array */
cooling_Redshifts =
(float *)mymalloc("cNz", cooling_N_Redshifts * sizeof(float));
}
bufsize = 0;
/* get size of the buffer */
MPI_Pack_size(cooling_N_Temp, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_nH, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_Temp, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_He, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_SolarAbundances, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_Redshifts, MPI_FLOAT, MPI_COMM_WORLD, &size);
bufsize += size;
MPI_Pack_size(cooling_N_Elements * EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD,
&size);
bufsize += size;
MPI_Pack_size(cooling_N_SolarAbundances * EL_NAME_LENGTH, MPI_CHAR,
MPI_COMM_WORLD, &size);
bufsize += size;
/* allocate memory for the buffer */
buffer = (char *)mymalloc("buffer", bufsize);
if (ThisTask == 0) {
position = 0;
MPI_Pack(cooling_Temp, cooling_N_Temp, MPI_FLOAT, buffer, bufsize,
&position, MPI_COMM_WORLD);
MPI_Pack(cooling_nH, cooling_N_nH, MPI_FLOAT, buffer, bufsize, &position,
MPI_COMM_WORLD);
MPI_Pack(cooling_Therm, cooling_N_Temp, MPI_FLOAT, buffer, bufsize,
&position, MPI_COMM_WORLD);
MPI_Pack(cooling_HeFrac, cooling_N_He, MPI_FLOAT, buffer, bufsize,
&position, MPI_COMM_WORLD);
MPI_Pack(cooling_SolarAbundances, cooling_N_SolarAbundances, MPI_FLOAT,
buffer, bufsize, &position, MPI_COMM_WORLD);
MPI_Pack(cooling_Redshifts, cooling_N_Redshifts, MPI_FLOAT, buffer, bufsize,
&position, MPI_COMM_WORLD);
for (i = 0; i < cooling_N_Elements; i++)
MPI_Pack(cooling_ElementNames[i], EL_NAME_LENGTH, MPI_CHAR, buffer,
bufsize, &position, MPI_COMM_WORLD);
for (i = 0; i < cooling_N_SolarAbundances; i++)
MPI_Pack(cooling_SolarAbundanceNames[i], EL_NAME_LENGTH, MPI_CHAR, buffer,
bufsize, &position, MPI_COMM_WORLD);
}
MPI_Bcast(buffer, bufsize, MPI_PACKED, 0, MPI_COMM_WORLD);
if (ThisTask != 0) {
position = 0;
MPI_Unpack(buffer, bufsize, &position, cooling_Temp, cooling_N_Temp,
MPI_FLOAT, MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, cooling_nH, cooling_N_nH, MPI_FLOAT,
MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, cooling_Therm, cooling_N_Temp,
MPI_FLOAT, MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, cooling_HeFrac, cooling_N_He,
MPI_FLOAT, MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, cooling_SolarAbundances,
cooling_N_SolarAbundances, MPI_FLOAT, MPI_COMM_WORLD);
MPI_Unpack(buffer, bufsize, &position, cooling_Redshifts,
cooling_N_Redshifts, MPI_FLOAT, MPI_COMM_WORLD);
for (i = 0; i < cooling_N_Elements; i++)
MPI_Unpack(buffer, bufsize, &position, cooling_ElementNames[i],
EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD);
for (i = 0; i < cooling_N_SolarAbundances; i++)
MPI_Unpack(buffer, bufsize, &position, cooling_SolarAbundanceNames[i],
EL_NAME_LENGTH, MPI_CHAR, MPI_COMM_WORLD);
}
/* free allocated memory */
myfree(buffer);
}
/*
* ----------------------------------------------------------------------
* This routine sets the value of solar metallicity using solar element
* abundances stored in the cooling tables.
* ----------------------------------------------------------------------
*/
void set_solar_metallicity(void) {
int i;
if (ThisTask == 0) printf("\n\nSolar abundances:\n\n");
/* set solar metallicity value */
for (i = 0, SolarMetallicity = 1; i < cooling_N_SolarAbundances; i++) {
if (ThisTask == 0)
printf("Solar abundance of %s = %g\n", cooling_SolarAbundanceNames[i],
cooling_SolarAbundances[i]);
if (strcmp(cooling_SolarAbundanceNames[i], "Hydrogen") == 0)
SolarMetallicity -= cooling_SolarAbundances[i];
if (strcmp(cooling_SolarAbundanceNames[i], "Helium") == 0)
SolarMetallicity -= cooling_SolarAbundances[i];
}
if (ThisTask == 0) printf("\nSolarMetallicity = %g\n", SolarMetallicity);
}
/*
* ----------------------------------------------------------------------
* Make pointers to element names
* ----------------------------------------------------------------------
*/
void MakeNamePointers() {
int i, j, sili_index = 0;
ElementNamePointers =
(int *)mymalloc("coolNele", cooling_N_Elements * sizeof(int));
SolarAbundanceNamePointers =
(int *)mymalloc("coolNele2", cooling_N_Elements * sizeof(int));
for (i = 0; i < EAGLE_NELEMENTS; i++) {
if (strcmp(ElementNames[i], "Silicon") == 0) sili_index = i;
}
for (i = 0; i < cooling_N_Elements; i++) {
SolarAbundanceNamePointers[i] = -999;
ElementNamePointers[i] = -999;
for (j = 0; j < cooling_N_SolarAbundances; j++) {
if (strcmp(cooling_ElementNames[i], cooling_SolarAbundanceNames[j]) == 0)
SolarAbundanceNamePointers[i] = j;
}
if (strcmp(cooling_ElementNames[i], "Sulphur") == 0 ||
strcmp(cooling_ElementNames[i], "Calcium") ==
0) /* These elements are tracked! */
ElementNamePointers[i] = -1 * sili_index;
else {
for (j = 0; j < EAGLE_NELEMENTS; j++) {
if (strcmp(cooling_ElementNames[i], ElementNames[j]) == 0)
ElementNamePointers[i] = j;
}
}
}
#ifdef EAGLE_VERBOSE
if (ThisTask == 0) {
for (i = 0; i < cooling_N_Elements; i++)
printf("cooling_ElementNames[%d] is cooling_SolarAbundanceNames %s\n", i,
cooling_SolarAbundanceNames[SolarAbundanceNamePointers[i]]);
printf("\n");
for (i = 0; i < cooling_N_Elements; i++)
printf("cooling_ElementNames[%d] is ElementNames %s\n", i,
ElementNames[ElementNamePointers[i]]);
printf("\n");
}
#endif
}
void eagle_set_abundances_to_solar() {
int i, element, index;
if (ThisTask == 0)
printf(" setting abundances of stars and particles to solar \n");
for (element = 0; element < cooling_N_Elements; element++) {
index = get_element_index(cooling_SolarAbundanceNames,
cooling_N_SolarAbundances, ElementNames[element]);
if (index >= 0) {
for (i = 0; i < NumPart; i++)
if (P[i].Type == 0) {
SphP[i].Metals[element] = cooling_SolarAbundances[element];
SphP[i].Metallicity = SolarMetallicity;
} else if (P[i].Type == 4)
StarP[P[i].pt.StarID].Metals[element] =
cooling_SolarAbundances[element];
}
}
}
/*
* ----------------------------------------------------------------------
* This routine initialize cooling at the beginning of the simulation
* ----------------------------------------------------------------------
*/
void eagle_init_cool() {
int i, j, k;
double Z_Solar_Calcium = 0, Z_Solar_Sulphur = 0, Z_Solar_Silicon = 0;
char fname[500];
/* use a general file to get the header */
if (ThisTask == 0) {
GetCoolingRedshifts();
sprintf(fname, "%sz_0.000.hdf5", All.CoolTablePath);
ReadCoolingHeader(fname);
}
//BroadCastHeader();
set_solar_metallicity();
if (element_present("Silicon") == 0) {
for (i = 0; i < cooling_N_SolarAbundances; i++) {
if (strcmp(cooling_SolarAbundanceNames[i], "Calcium") == 0)
Z_Solar_Calcium = cooling_SolarAbundances[i];
if (strcmp(cooling_SolarAbundanceNames[i], "Sulphur") == 0)
Z_Solar_Sulphur = cooling_SolarAbundances[i];
if (strcmp(cooling_SolarAbundanceNames[i], "Silicon") == 0)
Z_Solar_Silicon = cooling_SolarAbundances[i];
}
} else {
printf("Silicon is not present: Calcium and Sulphur will not be tracked!");
}
//cooling_ThermalToTemp =
// (float ****)mymalloc("coolTtoT", 2 * sizeof(float ***));
//for (i = 0; i < 2; i++) {
// cooling_ThermalToTemp[i] =
// (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
// for (j = 0; j < cooling_N_He; j++) {
// cooling_ThermalToTemp[i][j] =
// (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
// for (k = 0; k < cooling_N_nH; k++)
// cooling_ThermalToTemp[i][j][k] =
// (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
// }
//}
float cooling_ThermalToTemp[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
//cooling_SolarElectronAbundance =
// (float ***)mymalloc("coolSe", 2 * sizeof(float **));
//for (i = 0; i < 2; i++) {
// cooling_SolarElectronAbundance[i] =
// (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
// for (j = 0; j < cooling_N_nH; j++)
// cooling_SolarElectronAbundance[i][j] =
// (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
//}
float cooling_SolarElectronAbundance[2][cooling_N_nH][cooling_N_Temp];
/* Allocate for metal cooling */
//cooling_MetalsNetHeating =
// (float ****)mymalloc("CoolHeat", 2 * sizeof(float ***));
//for (i = 0; i < 2; i++) {
// cooling_MetalsNetHeating[i] =
// (float ***)mymalloc("coolHeat", cooling_N_Elements * sizeof(float **));
// for (j = 0; j < cooling_N_Elements; j++) {
// cooling_MetalsNetHeating[i][j] =
// (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
// for (k = 0; k < cooling_N_nH; k++)
// cooling_MetalsNetHeating[i][j][k] =
// (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
// }
//}
float cooling_MetalsNetHeating[2][cooling_N_Elements][cooling_N_nH][cooling_N_Temp];
/* Allocate for H and He cooling */
//cooling_HplusHeNetHeating =
// (float ****)mymalloc("CoolHHeHeat", 2 * sizeof(float ***));
//for (i = 0; i < 2; i++) {
// cooling_HplusHeNetHeating[i] =
// (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
// for (j = 0; j < cooling_N_He; j++) {
// cooling_HplusHeNetHeating[i][j] =
// (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
// for (k = 0; k < cooling_N_nH; k++)
// cooling_HplusHeNetHeating[i][j][k] =
// (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
// }
//}
float cooling_HplusHeNetHeating[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
//cooling_HplusHeElectronAbundance =
// (float ****)mymalloc("coolHHeNe", 2 * sizeof(float ***));
//for (i = 0; i < 2; i++) {
// cooling_HplusHeElectronAbundance[i] =
// (float ***)mymalloc("coolNHe", cooling_N_He * sizeof(float **));
// for (j = 0; j < cooling_N_He; j++) {
// cooling_HplusHeElectronAbundance[i][j] =
// (float **)mymalloc("coolNnH", cooling_N_nH * sizeof(float *));
// for (k = 0; k < cooling_N_nH; k++)
// cooling_HplusHeElectronAbundance[i][j][k] =
// (float *)mymalloc("coolNH", cooling_N_Temp * sizeof(float));
// }
//}
float cooling_HplusHeElectronAbundance[2][cooling_N_He][cooling_N_nH][cooling_N_Temp];
//MakeNamePointers();
/* Set H and He reionisation heating to erg/g */
All.REION_H_Heating_ERGpG = All.REION_H_Heating_EVpH * EV_TO_ERG / PROTONMASS;
All.REION_He_Heating_ERGpG =
All.REION_He_Heating_EVpH * EV_TO_ERG / PROTONMASS;
printf("Hydrogen reionisation heating = %e erg/g\n",
All.REION_H_Heating_ERGpG);
printf(" Helium reionisation heating = %e erg/g\n",
All.REION_He_Heating_ERGpG);
//cooling_ElementAbundance_SOLAR =
// (double *)mymalloc("coolNele", cooling_N_Elements * sizeof(double));
//cooling_ElementAbundance_SOLARM1 =
// (double *)mymalloc("coolNele", cooling_N_Elements * sizeof(double));
double cooling_ElementAbundance_SOLAR[cooling_N_Elements];
double cooling_ElementAbundance_SOLARM1[cooling_N_Elements];
#if EAGLE_TEST_COOL == 1
if (ThisTask == 0)
printf(
"EAGLE_TEST_COOL option 1: ThisTask = 0 will evaluate cooling rates, "
"then we will stop\n");
DriverEagleTestCool();
MPI_Barrier(MPI_COMM_WORLD);
printf("Test 1 of cooling finished\n");
endrun(-2);
#endif
#if EAGLE_TEST_COOL == 2
if (ThisTask == 0)
printf(
"EAGLE_TEST_COOL option 2: ThisTask = 0 will evaluate thermal "
"evolution at the mean density\n");
DriverEagleTestCoolThermalEvolution();
MPI_Barrier(MPI_COMM_WORLD);
printf("Test 2 of cooling finished\n");
endrun(-2);
#endif
}
/*
* ----------------------------------------------------------------------
* This routine allocates header arrays for the cooling tables
* ----------------------------------------------------------------------
*/
void allocate_header_arrays(void) {
int i;
/* bins for interpolation in log10(T [K]) */
cooling_Temp = (float *)mymalloc("coolTemp", cooling_N_Temp * sizeof(float));
/* bins for interpolation in hydrogen number density, log10(nh [/cm^3]) */
cooling_nH = (float *)mymalloc("coolNnH", cooling_N_nH * sizeof(float));
/* bins for interpolation in thermal energy per unit mass, log10(u [erg/g]) */
cooling_Therm = (float *)mymalloc("coolNT", cooling_N_Temp * sizeof(float));
/* bins for interpolation in Helium abundance */
cooling_HeFrac = (float *)mymalloc("coolNHe", cooling_N_He * sizeof(float));
/* table of solar abundances */
cooling_SolarAbundances =
(float *)mymalloc("coolNSol", cooling_N_SolarAbundances * sizeof(float));
/* assumed solar abundances used in constructing the tables, and corresponding
* names */
cooling_SolarAbundanceNames = (char **)mymalloc(
"coolNsolName", cooling_N_SolarAbundances * sizeof(char *));
if (ThisTask != 0)
for (i = 0; i < cooling_N_SolarAbundances; i++)
cooling_SolarAbundanceNames[i] =
(char *)mymalloc("coolNsolName2", EL_NAME_LENGTH * sizeof(char));
/* names of chemical elements stored in table */
cooling_ElementNames =
(char **)mymalloc("coolElNa", cooling_N_Elements * sizeof(char *));
if (ThisTask != 0)
for (i = 0; i < cooling_N_Elements; i++)
cooling_ElementNames[i] =
(char *)mymalloc("coolElNa2", EL_NAME_LENGTH * sizeof(char));
}
/*
* ======================================================================
* FUNCTIONS FOR BASIC TESTING OF COOLING ROUTINES - TEST_COOL flag.
* ======================================================================
*/
/*
* ----------------------------------------------------------------------
* This routine tests the interpolation of the cooling tables
* ----------------------------------------------------------------------
*/
#if EAGLE_TEST_COOL > 0
void DriverEagleTestCool() {
int iz;
int z_index;
float z1, z2, dz;
/* we use redshifts at exact values from the tables, */
/* as well as values in between */
for (iz = cooling_N_Redshifts - 2; iz >= 0; iz--) {
z1 = cooling_Redshifts[iz + 1]; /* collisional table */
z2 = cooling_Redshifts[iz]; /* collisional table */
get_redshift_index(z1, &z_index, &dz);
LoadCoolingTables(z_index);
/*
EagleTestCool(z1 - 0.1*(z1-z2));
*/
if (ThisTask == 0) EagleTestCool(z2);
}
}
void EagleTestCool(float z) {
FILE *outfile;
int i, j, k, index1, index2, H_index, He_index, z_index;
char filename[500];
double Hefrac, ener, dens, temp;
static double ParticleMetalFraction[EAGLE_NELEMENTS];
float dz;
/* given Helium abundance */
Hefrac = cooling_HeFrac[0];
get_redshift_index(z, &z_index, &dz);
H_index = element_index("Hydrogen");
He_index = element_index("Helium");
/* loop over elements */
for (i = 0; i < EAGLE_NELEMENTS; i++) {
/* put abundance of element i to solar value, zero all others */
for (j = 0; j < EAGLE_NELEMENTS; j++) ParticleMetalFraction[j] = 0;
index1 = get_element_index(cooling_ElementNames, cooling_N_Elements,
ElementNames[i]);
if (index1 >= 0) /* index < 0 for Helium which is not a metal */
{
index2 = get_element_index(cooling_SolarAbundanceNames,
cooling_N_SolarAbundances,
cooling_ElementNames[index1]);
if (index2 < 0 || index2 >= cooling_N_SolarAbundances) {
printf(" sorry: element %s not found \n", cooling_ElementNames[index1]);
endrun(-12344);
}
sprintf(filename, "testcool.z=%1.3f.Helium=%1.3f.Element=%s", z, Hefrac,
ElementNames[i]);
outfile = fopen(filename, "w");
fprintf(outfile, "# element = %s\n", ElementNames[i]);
ParticleMetalFraction[H_index] =
1.0 - Hefrac - cooling_SolarAbundances[index2];
ParticleMetalFraction[He_index] = Hefrac;
ParticleMetalFraction[i] = cooling_SolarAbundances[index2];
} else {
/* primordial cooling */
sprintf(filename, "testcool.z=%1.3f.Helium=%1.3f.Element=%s", z, Hefrac,
"Primordial");
outfile = fopen(filename, "w");
fprintf(outfile, "# element = %s\n", "Primordial");
ParticleMetalFraction[H_index] = 1.0 - Hefrac;
ParticleMetalFraction[He_index] = Hefrac;
}
/* write temperature density grid */
fprintf(outfile, "# redshift = %e\n", z);
fprintf(outfile, "# Helium abundance = %e\n", Hefrac);
fprintf(outfile, "# cooling_N_Temp = %d\n", cooling_N_Temp);
fprintf(outfile, "# cooling_N_nH = %d\n", cooling_N_nH);
fprintf(outfile, "# T[K] u[erg/g] N_h [cm^-3] Lambda [erg cm^3 s^-1]\n");
for (j = 0; j < cooling_N_Temp - 1; j++) {
ener = cooling_Therm[j];
for (k = 0; k < cooling_N_nH - 1; k++) {
dens = cooling_nH[k];
temp = eagle_convert_u_to_temp(dz, pow(10.0, ener), pow(10.0, dens),
Hefrac);
fprintf(outfile, " %e %e %e %e\n", temp, ener, pow(10.0, dens),
EagleCoolingRate(pow(10.0, ener), pow(10.0, dens), z, dz,
ParticleMetalFraction));
}
}
fclose(outfile);
}
}
/*
* ----------------------------------------------------------------------
* This routine tests the thermal evolution of gas at the critical density
* ----------------------------------------------------------------------
*/
void DriverEagleTestCoolThermalEvolution() {
FILE *outfile = 0;
char filename[100];
int i, HeIndex, HIndex;
int z_index;
float z, dz;
float zold, znow;
double zmax = 49.0, deltaz = 0.1, uinit = 5.0e10, delta = 0.0, rho, u;
double RhoCrit, hubble_z, dtime, da;
double element_metallicity[EAGLE_NELEMENTS];
double Hefrac = 0.248, Hfrac = 0.752, temp;
double uc = 0;
/* test only makes sense if co-moving integration is on */
if (All.ComovingIntegrationOn != 1) {
if (ThisTask == 0)
printf(
" this cooling test only works when All.ComovingIntegrationOn=1 \n");
endrun(-2);
}
for (i = 0; i < EAGLE_NELEMENTS; i++) element_metallicity[i] = 0;
HeIndex = element_index("Helium");
HIndex = element_index("Hydrogen");
element_metallicity[HeIndex] = Hefrac;
element_metallicity[HIndex] = Hfrac;
/* crit. density in code units */
RhoCrit = 3 * All.Hubble * All.Hubble / (8 * M_PI * All.G);
/* crit. density in physical g/cm^3 */
RhoCrit *= All.UnitDensity_in_cgs * All.HubbleParam * All.HubbleParam;
u = uinit;
if (ThisTask == 0) {
sprintf(filename, "thermevol.data");
if (!(outfile = fopen(filename, "w"))) {
printf("can't open file %s\n", filename);
endrun(888);
}
fprintf(outfile, "# delta = %g\n", delta);
fprintf(outfile, "# uinit = %g\n", uinit);
fprintf(outfile, "# metallicity = ");
for (i = 0; i < EAGLE_NELEMENTS; i++)
fprintf(outfile, " %g\t", element_metallicity[i]);
fprintf(outfile, "\n");
fprintf(outfile, "# z T[K] u[erg/g] rho[g/cm^3] reion heat [erg/g]\n");
}
for (z = zmax; z - deltaz > 0; z -= deltaz) {
znow = z - deltaz;
zold = z;
if (ThisTask == 0) {
printf(" REDSHIFT = %f, %f\n", zold, znow);
fflush(stdout);
}
u *= pow(((1 + znow) / (1 + zold)), 2);
get_redshift_index(znow, &z_index, &dz);
LoadCoolingTables(z_index);
/* physical baryon density at this z, for this over density, in g/cm^3 */
rho = (1 + delta) * RhoCrit * All.OmegaBaryon * pow(1 + znow, 3);
da = deltaz / pow(1 + znow, 2);
hubble_z = (All.HubbleParam * HUBBLE) *
sqrt(All.Omega0 * pow(1 + znow, 3) +
(1 - All.OmegaLambda - All.Omega0) * pow(1 + znow, 2) +
All.OmegaLambda);
dtime = (1 + znow) * da / hubble_z;
u = EagleDoCooling(u, rho, dtime, -deltaz, znow, dz, element_metallicity);
/* extra energy from Hydrogen reionization */
if (znow <= All.REION_H_ZCenter && znow + deltaz > All.REION_H_ZCenter) {
u += All.REION_H_Heating_ERGpG;
uc += All.REION_H_Heating_ERGpG;
if (ThisTask == 0)
printf(" Hydrogen reionzation: added %e erg/g at redshift %g\n", uc,
znow);
}
/* extra energy from Helium reionization */
/* Note that this term is added in the cooling calculation */
/* u += eagle_helium_reionization_extraheat(znow, -deltaz); */
uc += eagle_helium_reionization_extraheat(znow, -deltaz);
temp = eagle_convert_u_to_temp(dz, u, rho * Hfrac / PROTONMASS, Hefrac);
if (ThisTask == 0)
fprintf(outfile, " %g \t %g \t %g \t %g \t %g\n", z, temp, u, rho, uc);
}
if (ThisTask == 0) fclose(outfile);
}
#endif /* EAGLE_TEST_COOL > 0 */
#endif
src/cooling/EAGLE/eagle_interpol.h deleted 100644 → 0
+ 0
− 179
View file @ 261af9c3
#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 @ 89c22621
...@@ -68,15 +68,24 @@ int main() { ...@@ -68,15 +68,24 @@ int main() {
exit(1); 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_init(params, &us, &internal_const, &cosmo);
cosmology_print(&cosmo); cosmology_print(&cosmo);
cooling_init(params, &us, &internal_const, &cooling); cooling_init(params, &us, &internal_const, &cooling);
cooling_print(&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++){ for(int t_i = 0; t_i < n_t_i; t_i++){
...@@ -86,8 +95,6 @@ int main() { ...@@ -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.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)); p.entropy = pressure/(pow(p.rho,gamma));
gamma = 5.0/3.0;
cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const); cooling_du_dt = eagle_cooling_rate(&p,&cooling,&cosmo,&internal_const);
temperature_cgs = eagle_convert_u_to_temp(&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); 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