Skip to content
Snippets Groups Projects
Commit 52a289bd authored by Stefan Arridge's avatar Stefan Arridge
Browse files

Cooling works with single precision

parent d7a2da34
No related branches found
No related tags found
1 merge request!242Add the cooling infrastructure and the const_du and const_lambda
...@@ -3,30 +3,74 @@ import matplotlib.pyplot as plt ...@@ -3,30 +3,74 @@ import matplotlib.pyplot as plt
import h5py as h5 import h5py as h5
import sys import sys
filename = "./energy.txt" stats_filename = "./energy.txt"
plot_dir = "./" snap_filename = "uniformBox_000.hdf5"
cooling_lambda = 1.0e-22 #plot_dir = "./"
n_H = 0.1
#some constants in cgs units
k_b = 1.38E-16 #boltzmann
m_p = 1.67e-24 #proton mass
#initial conditions set in makeIC.py
rho = 3.2e6
P = 4.5e9
n_H_cgs = 0.1
gamma = 5./3.
T_init = 1.0e5 T_init = 1.0e5
T_floor = 1.0e4
plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H) #Read the units parameters from the snapshot
plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png" f = h5.File(snap_filename,'r')
#analytic_solution = np.zeros(n_snaps-1) units = f["InternalCodeUnits"]
array = np.genfromtxt(filename,skip_header = 1) 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)"]
parameters = f["Parameters"]
cooling_lambda = float(parameters.attrs["CreaseyCooling:Lambda"])
min_T = float(parameters.attrs["CreaseyCooling:minimum_temperature"])
mu = float(parameters.attrs["CreaseyCooling:mean_molecular_weight"])
X_H = float(parameters.attrs["CreaseyCooling:hydrogen_mass_abundance"])
#get number of particles
header = f["Header"]
n_particles = header.attrs["NumPart_ThisFile"][0]
#read energy and time arrays
array = np.genfromtxt(stats_filename,skip_header = 1)
time = array[:,0] time = array[:,0]
total_energy = array[:,2] total_energy = array[:,2]
total_mass = array[:,1]
#conversions to cgs
rho_cgs = rho * unit_mass / (unit_length)**3
time_cgs = time * unit_time
u_init_cgs = total_energy[0]/(total_mass[0]) * unit_length**2 / (unit_time)**2
#find the energy floor
u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
#find analytic solution
analytic_time = np.linspace(0.0,time_cgs[-1],10000)
du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
u_analytic = du_dt_cgs*analytic_time + u_init_cgs
cooling_time = u_init_cgs/(-du_dt_cgs)
#rescale energy to initial energy
total_energy /= total_energy[0]
u_analytic /= u_init_cgs
u_floor_cgs /= u_init_cgs
# plt.plot(e_kin_array,label = "Kinetic Energy") # plot_title = r"$\Lambda \, = \, %1.1g \mathrm{erg}\mathrm{cm^3}\mathrm{s^{-1}} \, \, T_{init} = %1.1g\mathrm{K} \, \, T_{floor} = %1.1g\mathrm{K} \, \, n_H = %1.1g\mathrm{cm^{-3}}$" %(cooling_lambda,T_init,T_floor,n_H)
# plt.plot(e_therm_array,label = "Internal Energy") # plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
plt.plot(time,total_energy,'r',label = "Total Energy") #analytic_solution = np.zeros(n_snaps-1)
#plt.plot(time,analytic_solution,'--k',label = "Analytic Solution") for i in range(u_analytic.size):
plt.xlabel("Time (code units)") if u_analytic[i]<u_floor_cgs:
plt.ylabel("Energy (code units)") u_analytic[i] = u_floor_cgs
#plt.ylim((0,1.0e8)) plt.plot(time_cgs,total_energy,'k',label = "Numerical solution")
plt.plot(analytic_time,u_analytic,'--r',lw = 2.0,label = "Analytic Solution")
plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
plt.plot((time_cgs[1],time_cgs[1]),(0,1),'m',label = "First output")
plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
plt.xlabel("Time (seconds)")
plt.ylabel("Energy/Initial energy")
plt.ylim(0,1)
plt.xlim(0,min(10.0*cooling_time,time_cgs[-1]))
plt.legend(loc = "upper right") plt.legend(loc = "upper right")
#plt.title(plot_title)
full_plot_filename = "%s/%s" %(plot_dir,plot_filename)
if (int(sys.argv[1])==0): if (int(sys.argv[1])==0):
plt.show() plt.show()
else: else:
......
...@@ -29,8 +29,8 @@ from numpy import * ...@@ -29,8 +29,8 @@ from numpy import *
periodic= 1 # 1 For periodic box periodic= 1 # 1 For periodic box
boxSize = 1 #1 kiloparsec boxSize = 1 #1 kiloparsec
L = int(sys.argv[1]) # Number of particles along one axis L = int(sys.argv[1]) # Number of particles along one axis
rho = 3.2e3 # Density in code units (0.1 hydrogen atoms per cm^3) rho = 3.2e6 # Density in code units (0.1 hydrogen atoms per cm^3)
P = 4.5e7 # Pressure in code units (at 10^5K) P = 4.5e9 # Pressure in code units (at 10^5K)
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 = "uniformBox.hdf5" fileName = "uniformBox.hdf5"
......
...@@ -56,12 +56,9 @@ void cooling_init(const struct swift_params* parameter_file, ...@@ -56,12 +56,9 @@ void cooling_init(const struct swift_params* parameter_file,
cooling->creasey_cooling.cooling_tstep_mult = parser_get_param_double(parameter_file, "CreaseyCooling:cooling_tstep_mult"); cooling->creasey_cooling.cooling_tstep_mult = parser_get_param_double(parameter_file, "CreaseyCooling:cooling_tstep_mult");
/*convert minimum temperature into minimum internal energy*/ /*convert minimum temperature into minimum internal energy*/
double k_b = phys_const->const_boltzmann_k; float u_floor = phys_const->const_boltzmann_k * cooling->creasey_cooling.min_temperature
double T_min = cooling->creasey_cooling.min_temperature; / (hydro_gamma_minus_one * cooling->creasey_cooling.mean_molecular_weight * phys_const->const_proton_mass);
double mu = cooling->creasey_cooling.mean_molecular_weight; float u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
double m_p = phys_const->const_proton_mass;
double u_floor = k_b * T_min / (hydro_gamma_minus_one * mu * m_p);
double u_floor_cgs = u_floor * units_cgs_conversion_factor(us,UNIT_CONV_ENERGY_PER_UNIT_MASS);
cooling->creasey_cooling.min_internal_energy = u_floor; cooling->creasey_cooling.min_internal_energy = u_floor;
cooling->creasey_cooling.min_internal_energy_cgs = u_floor_cgs; cooling->creasey_cooling.min_internal_energy_cgs = u_floor_cgs;
...@@ -95,7 +92,7 @@ void cooling_print(const struct cooling_data* cooling) { ...@@ -95,7 +92,7 @@ void cooling_print(const struct cooling_data* cooling) {
} }
void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us, void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us,
const struct cooling_data* cooling, struct part* p, double dt){ const struct cooling_data* cooling, struct part* p, float dt){
/*updates the entropy of a particle after integrating the cooling equation*/ /*updates the entropy of a particle after integrating the cooling equation*/
float u_old; float u_old;
...@@ -114,7 +111,7 @@ void update_entropy(const struct phys_const* const phys_const, const struct Unit ...@@ -114,7 +111,7 @@ void update_entropy(const struct phys_const* const phys_const, const struct Unit
/*This function integrates the cooling equation, given the initial /*This function integrates the cooling equation, given the initial
thermal energy, density and the timestep dt. Returns the final internal energy*/ thermal energy, density and the timestep dt. Returns the final internal energy*/
double calculate_new_thermal_energy(double u_old, double rho, double dt, float calculate_new_thermal_energy(float u_old, float rho, float dt,
const struct cooling_data* cooling, const struct cooling_data* cooling,
const struct phys_const* const phys_const, const struct phys_const* const phys_const,
const struct UnitSystem* us){ const struct UnitSystem* us){
......
...@@ -43,21 +43,21 @@ struct cooling_data { ...@@ -43,21 +43,21 @@ struct cooling_data {
#ifdef CONST_COOLING #ifdef CONST_COOLING
struct { struct {
double lambda; float lambda;
double min_energy; float min_energy;
double cooling_tstep_mult; float cooling_tstep_mult;
} const_cooling; } const_cooling;
#endif #endif
#ifdef CREASEY_COOLING #ifdef CREASEY_COOLING
struct { struct {
double lambda; float lambda;
double min_temperature; float min_temperature;
double hydrogen_mass_abundance; float hydrogen_mass_abundance;
double mean_molecular_weight; float mean_molecular_weight;
double min_internal_energy; float min_internal_energy;
double min_internal_energy_cgs; float min_internal_energy_cgs;
double cooling_tstep_mult; float cooling_tstep_mult;
} creasey_cooling; } creasey_cooling;
#endif #endif
}; };
...@@ -93,8 +93,8 @@ void cooling_init(const struct swift_params* parameter_file, ...@@ -93,8 +93,8 @@ void cooling_init(const struct swift_params* parameter_file,
void cooling_print(const struct cooling_data* cooling); void cooling_print(const struct cooling_data* cooling);
void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us, void update_entropy(const struct phys_const* const phys_const, const struct UnitSystem* us,
const struct cooling_data* cooling, struct part* p, double dt); const struct cooling_data* cooling, struct part* p, float dt);
double calculate_new_thermal_energy(double u_old, double rho, double dt, float calculate_new_thermal_energy(float u_old, float rho, float dt,
const struct cooling_data* cooling, const struct cooling_data* cooling,
const struct phys_const* const phys_const, const struct phys_const* const phys_const,
const struct UnitSystem* us); const struct UnitSystem* us);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment