diff --git a/examples/CoolingBox/coolingBox.yml b/examples/CoolingBox/coolingBox.yml
index 92cd2e733cc6f9e52ff750b3425d06d625732eeb..1fdebd74bfff9c187d99541f44ae0e5f2582185d 100644
--- a/examples/CoolingBox/coolingBox.yml
+++ b/examples/CoolingBox/coolingBox.yml
@@ -10,7 +10,7 @@ InternalUnitSystem:
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 4. # The end time of the simulation (in internal units).
- dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_min: 1e-4 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
@@ -36,7 +36,7 @@ InitialConditions:
# Dimensionless pre-factor for the time-step condition
LambdaCooling:
- lambda_cgs: 0.0 # Cooling rate (in cgs units)
+ lambda_cgs: 1.0e-22 # Cooling rate (in cgs units)
minimum_temperature: 1.0e4 # Minimal temperature (Kelvin)
mean_molecular_weight: 0.59 # Mean molecular weight
hydrogen_mass_abundance: 0.75 # Hydrogen mass abundance (dimensionless)
diff --git a/examples/CoolingBox/energy_plot.py b/examples/CoolingBox/energy_plot.py
index 57eb12d0dea44c8756a374b88c0188998d51f586..00e6fd1dfa0ee9bfbb9b5147282776f635b060f5 100644
--- a/examples/CoolingBox/energy_plot.py
+++ b/examples/CoolingBox/energy_plot.py
@@ -11,8 +11,8 @@ snap_filename = "coolingBox_000.hdf5"
k_b = 1.38E-16 #boltzmann
m_p = 1.67e-24 #proton mass
#initial conditions set in makeIC.py
-rho = 3.2e2
-P = 4.5e5
+rho = 3.2e3
+P = 4.5e6
#n_H_cgs = 0.0001
gamma = 5./3.
T_init = 1.0e5
@@ -32,51 +32,65 @@ X_H = float(parameters.attrs["LambdaCooling: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]
-total_energy = array[:,2]
+kin_plus_therm = array[:,2]
+radiated = array[:,6]
total_mass = array[:,1]
+#ignore first row where there are just zeros
time = time[1:]
-total_energy = total_energy[1:]
+kin_plus_therm = kin_plus_therm[1:]
+radiated = radiated[1:]
total_mass = total_mass[1:]
+total_energy = kin_plus_therm + radiated
+initial_energy = total_energy[0]
#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
+initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)**2
n_H_cgs = X_H * rho_cgs / m_p
+
#find the energy floor
u_floor_cgs = k_b * min_T / (mu * m_p * (gamma - 1.))
+
#find analytic solution
-analytic_time_cgs = np.linspace(time_cgs[0],time_cgs[-1],1000)
+analytic_time_cgs = np.linspace(0,time_cgs[-1],1000)
du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
-u_analytic = du_dt_cgs*(analytic_time_cgs - analytic_time_cgs[0]) + u_init_cgs
-cooling_time = u_init_cgs/(-du_dt_cgs)
+u_analytic_cgs = du_dt_cgs*analytic_time_cgs + initial_energy_cgs
+cooling_time_cgs = initial_energy_cgs/(-du_dt_cgs)
+
+for i in range(u_analytic_cgs.size):
+ if u_analytic_cgs[i]<u_floor_cgs:
+ u_analytic_cgs[i] = u_floor_cgs
+
+#rescale analytic solution
+u_analytic = u_analytic_cgs/initial_energy_cgs
#put time in units of cooling_time
-time=time_cgs/cooling_time
-analytic_time = analytic_time_cgs/cooling_time
-#rescale energy to initial energy
-total_energy /= total_energy[0]
-u_analytic /= u_init_cgs
-u_floor_cgs /= u_init_cgs
-# 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)
-# plot_filename = "energy_plot_creasey_no_cooling_T_init_1p0e5_n_H_0p1.png"
-#analytic_solution = np.zeros(n_snaps-1)
-for i in range(u_analytic.size):
- if u_analytic[i]<u_floor_cgs:
- u_analytic[i] = u_floor_cgs
-plt.plot(time-time[0],total_energy,'kd',label = "Numerical solution")
-plt.plot(analytic_time-analytic_time[0],u_analytic,'r',lw = 2.0,label = "Analytic Solution")
+time=time_cgs/cooling_time_cgs
+analytic_time = analytic_time_cgs/cooling_time_cgs
+
+#rescale (numerical) energy by initial energy
+radiated /= initial_energy
+kin_plus_therm /= initial_energy
+total_energy = kin_plus_therm + radiated
+plt.plot(time,kin_plus_therm,'kd',label = "Kinetic + thermal energy")
+plt.plot(time,radiated,'bo',label = "Radiated energy")
+plt.plot(time,total_energy,'g',label = "Total energy")
+plt.plot(analytic_time,u_analytic,'r',lw = 2.0,label = "Analytic Solution")
+#plt.plot(analytic_time,1-u_analytic,'k',lw = 2.0)
#plt.plot((cooling_time,cooling_time),(0,1),'b',label = "Cooling time")
#plt.plot((time[1]-time_cgs[0],time_cgs[1]-time_cgs[0]),(0,1),'m',label = "First output")
-plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
-plt.xlabel("Time (cooling time)")
-plt.ylabel("Energy/Initial energy")
-plt.ylim(0,1)
-plt.xlim(0,min(10,time[-1]))
+#plt.title(r"$n_H = %1.1e \, \mathrm{cm}^{-3}$" %n_H_cgs)
+plt.xlabel("Time / cooling time")
+plt.ylabel("Energy / Initial energy")
+#plt.ylim(0,1.1)
+plt.ylim(0.999,1.001)
+#plt.xlim(0,min(10,time[-1]))
plt.legend(loc = "upper right")
if (int(sys.argv[1])==0):
plt.show()
diff --git a/examples/CoolingBox/makeIC.py b/examples/CoolingBox/makeIC.py
index 79d7bd253999504471d0afc9bb1536e2829818ac..5de012a17af4eef71e56548602e7956faef529f5 100644
--- a/examples/CoolingBox/makeIC.py
+++ b/examples/CoolingBox/makeIC.py
@@ -29,8 +29,8 @@ from numpy import *
periodic= 1 # 1 For periodic box
boxSize = 1 # 1 kiloparsec
L = int(sys.argv[1]) # Number of particles along one axis
-rho = 3.2e6 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
-P = 4.5e9 # Pressure in code units (at 10^5K)
+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)
gamma = 5./3. # Gas adiabatic index
eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "coolingBox.hdf5"
diff --git a/examples/CoolingBox/run.sh b/examples/CoolingBox/run.sh
index 87a322465dcbfb439ee2d3f9de7df5ae109f8c72..c9a82469b3b03cf5ecc53aaa1b303f89ad40fc86 100755
--- a/examples/CoolingBox/run.sh
+++ b/examples/CoolingBox/run.sh
@@ -9,6 +9,6 @@ python makeIC.py 10
#-C 2>&1 | tee output.log
-#python energy_plot.py 0
+python energy_plot.py 0
-python test_energy_conservation.py 0
+#python test_energy_conservation.py 0
diff --git a/src/const.h b/src/const.h
index 1082bccc0bd514316c1457d67ea3f1023943c9e5..42c462bb836747393279ad384a9ee0193762f802 100644
--- a/src/const.h
+++ b/src/const.h
@@ -90,14 +90,14 @@
#define const_gravity_eta 0.025f
/* External gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
-//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
+//#define EXTERNAL_POTENTIAL_POINTMASS
+#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Cooling properties */
-#define COOLING_NONE
+//#define COOLING_NONE
//#define COOLING_CONST_DU
-//#define COOLING_CONST_LAMBDA
+#define COOLING_CONST_LAMBDA
//#define COOLING_GRACKLE
/* Are we debugging ? */
diff --git a/src/space.c b/src/space.c
index a9958f6fbd7d85060db99a9682b0de10f507085d..9cf534f8f58a3972fefa0db6042115a6edd4579f 100644
--- a/src/space.c
+++ b/src/space.c
@@ -1619,8 +1619,10 @@ void space_init(struct space *s, const struct swift_params *params,
} else {
for (size_t k = 0; k < Npart; k++)
for (int j = 0; j < 3; j++)
- if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j])
+ if (parts[k].x[j] < 0 || parts[k].x[j] >= dim[j]){
+ printf("parts[%lld].x[%d] = %f , dim[%d] = %f\n" , k , j , parts[k].x[j] , j , dim[j]);
error("Not all particles are within the specified domain.");
+ }
}
/* Same for the gparts */